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ABSTRACT 



The cores in molecular clouds are the densest and coldest regions of the interstellar medium (ISM). In these regions ISM-dust grains 
have the potential to coagulate. This study investigates the collisional evolution of the dust population by combining two models: a 
binary model that simulates the collision between two aggregates and a coagulation model that computes the dust size distribution 
with time. In the first, results from a parameter study quantify the outcome of the collision - sticking, fragmentation (shattering, 
breakage, and erosion) - and the effects on the internal structure of the particles in tabular format. These tables are then used as input 
for the dust evolution model, which is applied to an homogeneous and static cloud of temperature 10 K and gas densities between 10^ 
and 10' cm"'. The coagulation is followed locally on timescales of ~10' yr. We find that the growth can be divided into two stages: a 
growth dominated phase and a fragmentation dominated phase. Initially, the mass distribution is relatively narrow and shifts to larger 
sizes with time. At a certain point, dependent on the material properties of the grains as well as on the gas density, collision velocities 
will become sufficiently energetic to fragment particles, halting the growth and replenishing particles of lower mass. Eventually, a 
steady state is reached, where the mass distribution is characterized by a mass spectrum of approximately equal amount of mass per 
logarithmic size bin. The amount of growth that is achieved depends on the cloud's lifetime. If clouds exist on free-fall timescales the 
effects of coagulation on the dust size distribution are very minor. On the other hand, if clouds have long-term support mechanisms, 
the impact of coagulation is important, resulting in a significant decrease of the opacity on timescales longer than the initial collision 
timescale between big grains. 

Key words. ISM: dust, extinction - ISM: clouds - Turbulence - Methods: numerical 



1. Introduction 

Dust plays a key role in molecular clouds. Extinction of pen- 
etrating FUV photons by small dust grains allows molecules 
to survive. At the same time, gas will accrete on dust 
grains forming ice mantles consisting of si mple molecules 
jTielens & Hageniri982l; iHasegawa et al.lll992l) . Moreover, sur- 
face chemistry provides a driving force towards molecular 
complexity (Charnley et al. 1992; Aikawaetal. 2008). Finally, 
dust is often used as a proxy for the total gas (H2) col- 
umn density, either through near-IR extinc tion measurements 
or through sub-millimeter emission studies (Johnstone & Ballvl 



forming ice mantles. This growth is limited, however, because 
there are many small grains - which dominate the total grain 
surface area - over which the ice should be distributed; if all 
the condensible gas freezes out, the thickness of the ice mantles 
is still only 175 A (lDraineiri985h . Therefore, in dense clouds, 
coagulation is potentially a much more important promoter of 
dust growth. On a long timescale (>10'* yr), the interstellar grain 
size distribution is thought to reflect a balance between coagula- 
tion in dense clouds and shattering in interstellar shocks as ma- 
terial constantly cycles between dense and diffiise ISM phases 
(iJones et al.lll996l:(Dominik & Tielens|[T997h . 



l2006l;lAlves et al.l2OO7l:|j0rgensen et a l."2008'). Dust is often pre 
ferred as a tracer in these types of studies because it is now 
well established that - except for pure hydrides - all species 
condense out in th e form of ice mantles at the hig h densities 
of prestellar cores (Flower et al. ! 2006; B ergin & T afalla 2007; 
lAkvilmazet"ani2007r Thus, it is clear that our assessment of 
the molecular contents of clouds, as well as the overall state of 
the star and planet formation process, are tied to the properties 
of the dust grains - in particular, its size distribution. 

The properties of interstellar dust are, however, expected to 
change during its sojourn in the molecular cloud phase. First, 
condensation from the gas phase causes grain sizes to increase, 



Infrared and visual studies of the wavelength dependence 
of linear polarization and the ratio of total-to-selective extinc- 
tion were among the first observational indicati ons of the im- 
portance of grain growth in molecular clouds dCarrasco et al.l 
1973; Wilking et al. 1980; Whittet 2005). Early support for grain 
growth by coagulation in molecular clouds was also provided 
by a Copernicus study that revealed a decreased amount of vi- 
sual extinction per H-nucleus in the p-O ph cloud re lative to the 
value in the diffuse interstellar medium ('J uralll980l) . These type 
of visual and UV studies are by necessity limited to the outskirts 
of molecular clouds. Subsequent IR missions provided unique 
handles on the properties of dust deep inside dense clouds. In 
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particular, comparison of far-IR emission maps taken by IRAS 
and Spitzer and near-IR extinction maps derived from 2MASS 
indicate grain growth in the higher density regions ( Schn ee et alj 
I2OO8 ). Likewise, evidence for grain coagulation is also provided 
by a comparison of visual absorption studies (e.g., star counts) 
and sub-millimeter emission studies which imply that the small- 
est grains have been removed efficiently from the interstellar 
grain size distribution (Stepniket al. 2003). Similarly, a recent 
comparison of Spitzer-based, mid-IR extinction and submillime- 
ter emission studies of the dust characteristics in cloud cores re- 
veals systematic variations in the characteristics as a function of 
density co nsistent with models of grain growth by coagulation 
dButler & Ta n 2009). Dust-to-gas ratios derived from a compar- 
ison of line and continuum sub-millimet er data is also consisten t 
with grain growth in dense cloud cores dGoldsmith et al.|[l997l) . 
In recent years. X-ray absorption studies with Chandra have pro- 
vided a new handle on the total H column along a line of sight - 
that can potentially probe much deeper inside molecular clouds 
than UV studies - and in combination with Spitzer data, the de- 
creased dust extinction per H-nucleus reveals grain growth in 
molecular clouds (Winston et al. 2007). Finally, Spitzer/IRS al- 
lows studies of the silicate extinction inside dense clouds and 
a comparison of near-IR color excess with 10 jum optical depth 
reve als systematic vari ations which is likely caused by coagula- 
tion dChiar et aLll2007h . This is supported by an analysis of the 
detailed absorption profi le of the 10 silicate absorption band 
in the se environments dBowev et al.l Il998l: Ivan Breemen et alj 
l2009h . 

Because it is the site of planet formation, theoretical coag- 
ulation studies have largel y focused on grain gr owth in pro- 
toplanetary disks ( Weide nschilUng & Cuzziiri993h . In molecu- 
lar clouds, dust co agul ation has been theoretically modeled b y 
lOssenkopj dl993h and IWeidenschilling & Ruzmaikinal dl994l) . 
In these studies, coagulation is driven by processes that pro- 
vide grains with a relative motion. For larger grains (>100 A) 
turbulence in particularly is important in providing relative ve- 
locities. These motions - and the outcomes of the collisions - 
are very sensitive to the coupling of the particles to the turbu- 
lent eddies, which is determined by the surface area-to-mass ra- 
tio of the dust particles. At low velocities, grain collisions will 
lead to the growth of very open and fluffy structures, while at in- 
termediate velocities compaction of aggregates occurs. At very 
high velocities, cratering and catastrophic destruction will halt 
the growth CDominik & Tielens 1997; Paszun & Dominik 2009; 
iBlum & Wurml2008h . Thus, to study grain growth requires us to 
understand the relationships between the macroscopic velocity 
field of the molecular cloud, the internal structure of aggregates 
(which follows from its collision history), and the microphysics 
of dust aggregates collisions. In view of the complexity of the 
coagulation process and the then existing, limited understand- 
ing of the coagulation process, previous studies of coagulation 
in molecular cloud settings have been forced to make a num- 
ber of simplifying assumptions concerning the characteristics of 
growing aggregates. 

Theoretically, our understanding of the coagulation process 
has been much improved by the development of the atomic force 
microscope, which has provided much insight in the binding of 
individual monomers. This has been translated into simple re- 
lationships between velocities and material parameters, which 
prescribe under whic h conditions sticking, compaction, and 
fragmentation occur (^ Chokshi et alj ll993l ; lDominik & TielensI 
[l997). Over the last decade, a number of elegant experimen- 
tal studies by Blum and coworkers (see iBlum & WurmI |2008|) 
have provided direct support for these concepts and in many 



ways expanded our understanding of the coagulation process. 
Numerical simulations have translated these concepts into sim- 
ple recipes, which link the collisional parameters and the ag- 
grega te properties to the stru ctures of the evolving aggre- 
gates dPaszun & Dominikl2009l) . Together with the development 
of Monte Carlo methods, in which particles are indiv idually 
followed dOrmel et all 120071: IZsom & Dull emond"2008). these 
studies provide a much better framework for modeling the coag- 
ulation process than hitherto possible. 

In this paper, we reexamine the coagulation of dust grains in 
molecular cloud cores in the light of this improved understand- 
ing of the basic physics of coagulation with a two-fold goal. 
First, we will investigate the interrelationship between the de- 
tailed prescriptions of the coagulation recipe and the structure, 
size, and mass of aggregates that result from the collisional evo- 
lution. Therefore, a main goal of this work is to explore the full 
potential of the collision recipes, e.g., by running simulations 
that last long enough for fragmentation phenomena to become 
important. Second, we aim to give a simple prescriptions for 
the temporal evolution of the total grain surface area in molecu- 
lar clouds, thereby capturing its observational characteristics, in 
terms of the physical conditions in the core. 

This paper is organized as follows. Section|2]presents a sim- 
plified and static model of molecular clouds we adopted in our 
calculations. Section[3]describes the model to treat collisions be- 
tween dust grains and, more generally, aggregates of dust grains. 
In Sect. |4] the results are presented: we discuss the imprints of 
the collision recipe on the coagulation and also present a param- 
eter study, varying the cloud gas densities and the dust material 
properties. In Sect. |5] we review the implications of our result to 
molecular clouds and Sect. |6] summarizes the main conclusions. 



2. Density and velocity structure of molecular 
clouds 

The physical structure of molecular clouds - the gas density and 
temperature profiles - is determined by its support mechanisms: 
thermal, rotation, magnetic fields, or turbulence. If there is only 
thermal support to balance the cloud's self-gravity and the tem- 
perature is constant, the density, assuming spherical symmetry, 
is given by pg oc r"^, where pg is the gas density and r the ra- 
dius0 However, the isotherma l sphere is unstable as it heralds 
the collapse phase dShulll977h . The cloud then collapses on a 
free-fall timescale 



%=,^ = l.lxlO=yrf— 

^32Gpg ^ \105cm-3j 



(1) 



where G is Newton's gravitational constant, n - pg/mufJ^ the 
number density of the molecular gas|3 nin the hydrogen mass, 
and fi - 2.34 the mean molecular mass. Thermally supported 
cores are stable if the thermal pressure wins over gravity, a sit- 
uation described by the Bonnor-Ebert sphere, where an external 
pressure confines the cloud (still assuming a constant tempera- 
ture). 

Magnetic fields in particular are important to support the 
cloud against the opposing influence of gravity. Ion-molecule 
collisions provide friction between the ions and neutrals and in 



' A list of symbols is provided in Appendix IdI 

- Note that our definition of density - n, the number of gas molecules 
per cm' - differs from the density of hydrogen nuclei (iin), which is 
more commonly referred to as density in the dust/ISM communities. 
For cosmic abundances nn is related to n as ;i ^ 1.7nH- 
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that way couple the magnetic field to the neutral cloud. Over 
time the magnetic field will slowly leak out on an ambipolar dif- 
fusion timescale 

(«)~^''^^^° ^'^(lO^cm-O ' 

where Ki„ is the ion-molecular collision rate and n, the density 
in ions. In Eq. (|2]) we have used an ion-neutral collision rate of 
Ki„ = 2 X 10"^ cm^ s ' and a degree of ion ization due to cosmic 
rays of n,/n = 2 x lO^VV" (lTielensll2005"l) . 

Turbulence is another possible support mechanism of molec- 
ular cores, but its nature is dynamic - rather than (quasi)static. 
At large scales it provides global support to molecular clouds, 
whereas at small scales it locally compresses the gas. If these 
overdensities exist on timescales of Eq. ([T]i, collapse will follow. 
This is the gravo-turbulent fragmentation picture of turbulence- 
dominated molecular clouds, where the (supersonic) turbulence 
is driven at large scal es, but also reaches the scal es of quies- 
cent ( subsonic) cores ( Mac Low & Klessen 2004; Klessen et alj 
l2005h . In this dynamical, turbulent-driven picture both molecu- 
lar clouds and cores are transient objects. 

Thus, cloud cores will dynamically evolve due to either am- 
bipolar diffusion and loss of supporting magnetic fields or due to 
turbulent dissipation, or simply because the core is only a tran- 
sient phase in a turbulent velocity field. In this work, for reasons 
of simplicity, we consider only a static cloud model - the work- 
ing model - in which turbulence is unimportant for the support 
of the core, but where (subsonic) turbulence is included in the 
formalism for the calculation of relative motions between the 
dust particles. 

2.1. Working model 

In this exploratory study we will for simplicity adopt an homo- 
geneous core of mass given by the critical Jeans mass. Moreover, 
we assume the cloud is turbulent, but neglect the influence 
of the turbulence on the support of the cloud. Thus, our ap- 
proximation is probably applicable for high density, low mass 
cores as veloci ty dispersions increase towards high mass cores 
(iKawamura et al.l 119981) . The homogeneous structure ensures 
that collision timescales are the same throughout the cloud, i.e., 
the coagulation and fragmentation can be treated locally. In our 
calculations, the sensitivity of the coagulation on the gas density 
n will be investigated and the relevant coagulation and fragmen- 
tation timescales will be compared to the timescales in Eqs. ([B 
and ©. 

Starting from the isod ense sphere, the cloud out er radius is 
given by the Jeans length (iBinnev & Tremain3ll987l) 

1 Ucl I n I T 

^-Hci^o-^^^p'^d^) (loi) ' 

where Cg = kT 1^ is the isothermal sound speed and T the tem- 
perature of the cloud. A temperature of 10 K is adopted. The 
sound crossing time Lj/cg is comparable to the free-fall time of 
the cloud. 

Regarding the driving scales of the turbulence we assume 
(i) that the largest eddies decay on the sound crossing time, i.e., 
fL = Lj/cg, and (ii) that the fluctuating velocity at the largest 
scale is given by the sound speed, = Cg. Thus, the turbulent 
viscosity is Vt - Lu]^ = u^ti = CgLj with L = Lj the size of the 
largest eddies. Although our parametrization of the large eddy 



quantities seems rather ad-hoc, we can build some trust in this 
relation by considering the energy dissipation rate v^/L per unit 
mass, which translates into a heating rate of 

nr = ipg . 2.5 X 10- erg cm- s- (^) (j^f ■ 

(4) 

Based upon observational studies of turbulence in cores. iTielensI 
(2005) gives a heating rate of nF = SxlO^^'^nerg s"', with which 
Eq. (IHi reasonably agrees for the range of densities we consider. 
Additionally, the adoption of the crossing time and sound speed 
for the large eddy properties are natural upper limits. 

The turbulent properties further follow from the Reynolds 
number. 

Re = A = = 6.2X10^ {^^ri^r, (5) 

Vm CgCfp/3 \105cm-3/ llOK/ 

where Vm is the molecular viscosity and f^fp the mean free path 
of a gas particle. Assuming a Kolmogorov cascade, the turn-over 
time and velocity at the inner scale follow from the Reynolds 
number 

/ „ \-3/4 / ■r -,-1/4 

, = Re-/V.2.2xlO'yr(^^^) (— ) ; (6a) 



V, = Re-^'W = 2.1 X 10^ cm s'' — . 

llO^cms-i/ llOK/ 

(6b) 

A collisional evolution model requires a prescription for the 
relative velocities Ai> between two solid particles. Apart from 
turbulence, other mechanisms, reflecting differences in the ther- 
mal, electrostatic, and aerodynamic properties of particles, will 
also provide particles with a relative motion. However, under 
most molecul ar cloud condition s turbulence will dominate the 
velocity field dOssenkopSl 19931) and in this work we only con- 
sider turbulence. Then, the surface area-to-mass ratio of the dust 
particles is of critical importance since this quantity determines 
the amount of coupling between the du st particles and the gas . 
We use the analytic approximations of lOrmel & Cuzzil (|2007|) 
for the relative velocity between two particles. These expres- 
sions only include contributions that arise as a result of the parti- 
cle's inertia in a turbulent velocity field and do not contain con- 
tributions that a rise from gyroresonance acceleration (see, e.g., 
lYan et al.ll2004l) . See Appendix|A]for more details. 

3. Collision model 

Dust grains that collide can stick together, forming aggregates 
(see Fig. [T]). In this work we consider the collisional evolution of 
the distribution of aggregates, f(N, t): the number of aggregates 
of mass with time. Many works have studied aggregate growth 
under the conditions of perfect sticking upon contact, neglect- 
ing th e effects of the impact energy on the structure of aggre - 
gates (lMeakinlll988l: iMeakin & Donnlll988t lbssenkopj|I993h . 
Of special interest are the particle-cluster aggregation (PCA) and 
cluster-cluster aggregation (CCA) modes. In PCA, the aggre- 
gates collide only with single grains, while in CCA the collision 
partners are of similar size and structure. In CCA, the emerging 
structures become true fractals, with a fractal dimension ~2. For 
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Fig. 1. 2D projection of a fluffy aggregate with indication of the 
geometrical radii, Oq-, and the outer radii, flout- 

PCA, however, an homogeneous structure is eventually reached 
at a filling factor of ~15%. 

However, the assumption that the internal structure is fixed 
(as in fractals) becomes invalid if the collisions take place be- 
tween particles of diff'erent size. Furthermore, at higher energies 
the assumption of 'hit-and-stick' breaks down: aggregate bounc- 
ing, compaction (in which the constituent grains rearrange them- 
selves), and fragmentation lead to a rearrangement of the inter- 
nal structure. These collisional processes, except bouncing, are 
included in our collision model. 

We quantify the internal structure of aggregates in terms of 
the geometrical filling factor, <pa-, defined as 




where we have assumed that the aggregate contains equal size 
grains of radius flo with flo- = y/cr/n the projected surface equiv- 
alent radius of the aggregate. For very fluffy aggregates flo- can be 
much less than the outer radius of the aggregate, flout, see Fig.[T] 
The definition of the filling factor in terms of the projected area 
determines its aerodynamic behavior, and thereby the relative 
velocities (An) between the dust aggregates^ 

Each collisions is classified into one of three groups: 

1. Hit-and-stick. At low collision energies, the internal struc- 
ture of the particles is preserved. 

2. Local. Only a small part of the aggregate is affected by the 
collision, as in, e.g., erosion. The mass ratio between the two 
particles is large. 

3. Global. The collision outcome results in a major change to 
the structure or size of the target aggregate. Relevant for 
equal-size particles or at large energies. 

Figured provides the adopted decision chain between the three 
regimes. Parameters that enter the chain are the collision energy, 

E^U(A.f^'-^!^(A.f, (8) 
2 2 nil + m2 

^ The compactness para meter 0n- is the inv e rse of the enlarge ment 
factor if/, previously used in lOrmel et al-l ( l2007l) - [Ossenkopil ( Il993h uses 
X = (p'^^ as its structural parameter 



where the reduced mass, the particle masses (mi and m2, or, 
in number of grains, .^Vi and N2), and the critical energies Ei-ow 
and /ibr- Here, E\,i- is the energy to break a bond between two 
grains (the contact) and E^on the energy required to roll the con- 
tact area over a visi ble fraction of the grain- T hese critical ener- 
gies are defined as jDominik & Tieleiisll 19971) 

£br=Ato 57^; (9a) 

2*2/3 

E,oll = bTT^^cntyOfi, (9b) 

where fl^ is the reduced radius of the grains (fl^ = flo/2 for equal 
size monomers), 7 the surface energy density of the material, 
and £* the reduced elastic modulus- Here, following laboratory 
experiments (Blum & Wurm 2000J we adopt the values ^crft = 
2 X 10"^ cm and Abr = 2.8 x 10^, which are larger th an the 
theoretically derived values of iDominik & TielensI (Il997[) . 

Thus, when collisional energies are low enough for aggre- 
gate restructuring to be unimportant (experimentally determined 
to be Siiroii; B lum & Wurm 2000) particles are in the hit-and- 
stick regime. Similarly, when the collision energy is sufficient to 
break all contacts the collision falls - obviously - in the global 
regime. In the remainder the mass ratio of the colliding particles 
determines whether the collision is global or local. For mass ra- 
tios smaller than 0.1 the changes become more localized and it 
is seen from the simulations that at this point the energy dis- 
tribution during a collision becomes inhomogeneous. Thus, we 
take N2/N1 = 0-1 as the transition point- A further motivation for 
adopting this mass ratio is that we construct the global recipe out 
of simulations between aggregate of the same size- Therefore, 
the mass range which it represents should not become too large. 

Although in our collision model aggregates are character- 
ized by only two properties (A^ and (pa-), the collision outcome 
involves many other parameters (discussed in Sect. 13.31 ). These 
parameters are provided in tabulated form as a function of 
three input parameters - dimensionless energy parameter, fill- 
ing factor, and impact parameter b - for both the local and the 
global recipe. To obtain these collision parameters, direct nu- 
merical simulations between two colliding aggregates were per- 
formed, in which the equations of motions for all gr ains within 
the two colliding aggregates are computed jPaszun & DominiJ 
2009). An example of these quantities is the fraction of miss- 
ing collisions, which is a result of the fact that we have de- 
fined the collision cross section cr'- in terms of the outer radii, 
cr'^ - 7r(flout,i +flout,2)^- Appendix iBjpresents a description of the 
numerical simulations with their results, discusses a few auxil- 
iary relations that are required for a consistent treatment of the 
collision model, and treats the format of the collision tablesQ 

Two key limitation of these binary aggregate simulations fol- 
low from computational constraints: (i) the number of grains that 
can be included is limited to ~ 10^; and (ii) the simulations 
cannot take into account a grain size distribution that spans over 
orders of magnitude in mass. These limitations are reflected in 
our collision model and constitute a potential bottleneck for the 
level of realism for the application of our results to molecular 
clouds. We therefore first motivate our choice for the monomer 
size and present scaling relations that provide a way to extrapo- 
late the results beyond the parameter space sampled in the sim- 
ulation. 



* The tables are provided online. 
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Dust properties 

material properties y, £* 
grain size a^ 



Velocity field 



Particle #1 

mass Nj; filling factor ^ 



Particle #2 

mass filling factor j 



Critical energies 

breal<ing energy Ebr - A^, 



(£*)2/3 

rolling energy E,„n - 3;r^^eritr«o 



Collision properties 

collision energy E = -mnNui^vf 
2 ^ 

mass ratio l^/Ni 




Fig. 2. Schematic decision chain employed to distinguish between the hit-and-stick, global, and local recipes. 



3.1. Representative monomer size of the MRN grain 
distribution 

Our recipe is based on simulations of aggregates that are built of 
monomers of a single size. Therefore, we treat a monodisperse 
distribution of grains. In reality, interstellar dust exhibits a size 
distribution, or a seri es of size distribut i ons based on the vari- 
ous grain types (e.g.. iDesert et al.lll990l : IWeingartn er & Drairi^ 
I2OOI ). For simplicity, we compare our monodisperse approach 
with the MRN grain distribution, in which the number of grains 
decreases as a -7/2 power-law of size, i.e., n{a)da cx g'^^^, be - 
tween a lower (oi) and an upper (of) size (iMathis etal.lfT977l) . 
Thus, in the MRN-distribution the smallest grains dominate by 
number, whereas the larger grains dominate the mass. For an 
MRN distribution we adopt, Oi = 50 A and Of - 0.25 fim. To 
answer the question which grain size best represents the MRN 
distribution, we consider both its mechanical and aerodynamic 
properties. 

In the monodisperse situation the mechanical properties of 
a particle (its strength) can be estimated from the total binding 
energy per unit mass, Ebi /mo, if we assume each grain has one 
unique contact. In the literature the strength of a material is usu- 
ally denoted by the quantity Q. Thus, for a monodisperse model 
we have 



Go 



£br(flo/2) , (00/2)4/3 



k2 



"O ' 



(10) 



where k = 3A\„-'y^^^ /Anp^S*^^^ is a material constant with the 
bulk density of the material. A smaller grain size thus lead to 
significantly stronger aggregates. For the MRN distribution we 
assume that a typical contact always involves a small grain, i.e., 
G/j ^ a, enters in the E\,r expression. Assuming again that the 
number of contacts is of the order of the number of grains, their 
average strength is given by 



MRN - 



Ebi (oi) n(a)da 
4;rps/3 J"^ n(a)a^da 



(11) 



where we have used that Uf » a\. Equating Eqs. ( fTOt and ( fTTT i 
it follows that the grain size at which the monodisperse model 
gives aggregates that have the same strength as the MRN is oq - 
24/5flJ/"V/'° = 560A. 



Apart from the mechanical properties, the aerodynamic 
properties of aggregates are of crucial importance to the col- 
lisional evolution. This mainly concerns the initial (fractal) 
growth stage. For a single grain cr/m - (3/4ps)ag'. For the 
MRN distribution an upper limit on cr/m is provided by assum- 
ing that all of its surface is exposed, i.e., as in a 2D arrangement 
of grains; then. 



cr 
m 



Jj' 7Tn{a)a^da 3 
£' iAnpjy)n{a)aMa ~ M 



■(aiaf) 



-1/2 



(12) 



and the equivalent aerodynamic grain size becomes yfoiOf = 350 
A. However, this 2D result for the equivalent monodisperse size 
of the MRN distribution is a considerable underestimation, for 
three reasons: (i) in 3D the grains will overlap and cr becomes 
lower at the same mass; (ii) due to their low rolling energies, the 
smallest grains of size a, will already initiate restructure at very 
low velocities; (iii) in the case of ice-coating, the lower grain 
size fl, will be larger by a factor of ~4. 

In three dimensions, however, the definition of an equiva- 
lent aerodynamic size becomes ambiguous, because cr/m is not 
a constant. To nevertheless get a feeling of the trend, we have 
calculated the aerodynamic properties of MRN aggregates that 
consists out of a few big grains, such that their total compact 
volume is equivalent to a sphere of 0.2 - 0.3 /im. These MRN 
aggregates were constructed through a PCA process, i.e., adding 
one grain at a time. Because the aggregates contains the large 
grains, they fully sample the MRN distribution and can therefore 
be regarded as the smallest building blocks for the subsequent 
collisional evolution. 

We observed that, due to the above mentioned self-shielding, 
the aerodynamic size increases to ~0.08 - 0.1 2 yum, si gnificantly 
higher than the 2D limit of Eq. (fT2] i (see also lOssenkopf . 19931) . 
Thus, the initial clustering phase of MRN-grains produces struc- 
tures that behave aerodynamically as compact grains of ~0. 1 yum. 
We remark that this estimate is approximate - a CCA-like clus- 
tering will decrease it, whereas the above mentioned preferential 
compaction of the very small grains will increase cr/m - but the 
trend indicates that in 3D the aerodynamic size becomes skewed 
to the larger grains in the distribution. Therefore, we take 0. 1 fim 
as the equivalent monomer grain size of the MRN-distribution, 
but to assess the sensitivity of the adopted grain size to the results 
we also consider models with a different grain sizes. 



6 



C.W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. 



3.2. Scaling of the results 

A key limitation of the aggregate-aggregate collision simulations 
is the number of grains that can be used; typically, < 10^ is 
required in order to complete a thorough parameter study within 
a reasonable timeframe. As a consequence, the mass ratio of the 
colliding aggregate is also restricted. Furthermore, in the numer- 
ical experiments all simulations were performed using material 
properties applicable to silicates, whereas in molecular clouds 
we expect the grains to be coated with ice mantles. Clearly, scal- 
ing of the results is required such that the findings of the numer- 
ical experiments can be applied to aggregates of different size 
and composition. 

Therefore, we scale the collisional energy E to the critical 
energies, fibr and Eroii, since these quantities involve the mate- 
rial properties. For a collision between silicate aggregates and 
ice-coated aggregates a similar fragmentation behavior may be 
expected if the collisional energy in the latter case is a factor 
E!^^^/E^^_ higher Similarly, restructuring is determined by the 
rolling energy, Ei-on. Thus, the collision energy is normalized to 
£^roii where it concerns the change in filling factor and to Eb^ for 
all other parameters that quantify the collision outcome. 

The division between the global and local recipes is also 
closely linked to scaling arguments. In the global recipe ener- 
gies are normalized to the total number of monomers, A^tot- Thus, 
a collision taking place at twice the energy and twice the mass 
leads to the same fragmentation behavior However, in the local 
recipe the amount of damage will be independent of the size of 
the bigger particle. In this case we then scale by A'^, essentially 
the mass of the projectile. This information is captured in a sin- 
gle dimensionless energy parameter s. 



where E^h and A^eff depend on the context: the energy Eciit can 
be either one of E^^ or iiroii, whereas A^eff is one of A^tot or A'^ (see 
Tablel331l. 

3.3. The collision parameters 

In discussing the collision outcomes, we focus on the local and 
global recipes, which are a direct result of the numerical simu- 
lations. The hit-and-stick recipe is discussed in Appendix lB.2.31 
To streamline the recipe for a Monte Carlo approach, the speci- 
fication of the collision outcomes slightly deviates from our pre- 
vious study (Paszun & Dominik 2009). 

In the general case of a collision including fragmentation the 
emergent mass distribution, /(m), consists of two components: 
(i) a power-law component that describes the small fragments; 
and (ii) a large fragment component that consist out of one or 
two fragments (see Fig. [3]). The separation between the two com- 
ponents is set, somewhat arbitrarily, at a quarter of the total mass 
of the aggregates, A^tot - Ni + N2. (It turns out that for our simu- 
lations the precise place of the cut is unimportant, because of the 
lack of severe fragmentation events). Then, the power-law dis- 
tribution spans the range from monomer mass up to A^ = Mot/4, 
whereas the large-fragment component consists of zero, one, or 
two aggregates of masses larger than Mtot/4. To obtain the num- 
ber of large fragments, the recipes provide the mean number of 
large fragments, Nf, together with its spread Sf. 

Table [33] lists all quantities describing a collision outcome. 
Apart from Nf and Sf these include: 

- The fraction of missing collisions, /miss- This number gives 
the fraction of collisions in which no interaction between the 




1 



N 

Fig. 3. Sketch of the adopted formalism for the size distribution 
with which the results of the aggregate collision simulation are 
quantified. See text and Table 13. 3] for the description of the sym- 
bols. If /p„i = no power-law component exist. The Nf and Sf 
parameters essentially indicate whether we have zero, one, or 
two large fragments. 



aggregates took place. Missing collision are a result from the 
choice of normalizing the impact parameter b to the outer 
radius Oout, b - b/(aout,i + Oout.a) (see Appendix lB.2.2l i. For 
large values of b and very fluffy structures /Iniss becomes sig- 
nificant. 

- The mass fraction in the power-law component, /pwi. It gives 
the fraction of the total mass (Mot) that is contained in the 
power-law component. In the local recipe /pwi is defined rel- 
ative to M^,, because here the amount of erosion is measured 
with respect to the smaller projectile. 

- The exponent of the power-law distribution, q. It determines 
the distribution of the small fragments, i.e., /(m) oc m^''. 

- The relative change in filling factor, C^. It gives the change in 
filling factor of the large fragment component, (^"'^* - C^^™'. 

< I reflects compaction, whereas > 1 reflects decom- 
paction. Because refers to the chance in the filling factor 
of the large aggregate (for both the local and global recipe), 
its dimensionless energy parameter e is always normalized 
to the total number of grains, Meff = Mtot. Thus, the com- 
paction may be local and moderate, but the affected quantity 
- the filling factor - describes a global property. Moreover, to 
prevent possible spuriously high values of 0o-, we artificially 
assign an upp er limit of 33% to the co llisional compaction 
of aggregates dBlum & Schraplerfl2004l) . 

3.3.1. The local recipe 

Figure]?^ shows how much mass is ejected during collisions 
at different energies and for different filling factors (symbols). 
Recall that in the local recipe the /p„i quantity involves a nor- 
malization to N/j, rather than Mtot- At high energies, then, the 
excavated mass may exceed the mass of the small projectile by 
even two orders of magnitude. The distribution of the small frag- 
ments created by the erosion is relatively flat with slopes oscil- 
lating between q = -2.0 and q = -1.3. The number of large 
fragments Np rarely increases above unity. The exception are the 
'lucky projectiles' that destroy the central contacts of very fluffy 
aggregates, causing the two sides of the aggregate to become dis- 
connected. If energies are sufficiently high, fragments produced 
in a cratering event can result in secondary impacts, enhancing 
the erosion efficiency. 
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Symbol Description Energy scaling parameters in Eg. (I13"F ~ 

Global Local 

0) (2) (3) (4) 

/miss Fraction of collisions that resulted in a miss" — — 

A'f Mean number of large fragments Mot^^br N^E^r 

Sf Standard deviation of the A^f Mot^^br N^E-^^ 

/pwi The fraction of the mass in the small fragments component. Mot^^br N^E-^r 

Normalized to A^tot (global recipe) or A'^ (local recipe). 

q Exponent of the power-law distribution of small fragments Mot^^br N^E\,^ 

Cij, = ^a-l<t'™^ Relative change of the geometrical filling factor Mot^roU Mot^^roii 



Note. Columns (3)-(4) denote the energy scaling expressions used to obtain the dimensionless energy parameter, e, see Eq. l |13t . " Given as 
function of a„ut/acr instead of (pa-, see Appendix lB.3l 



Since the influence of the impact is local, the change in filling 
factor is relatively minor (see Fig. lib). However, increasing the 
collision energy results in an increasing degree of compression. 
Only very fluffy and elongated aggregates may break in half, 
causing an artificial increase of the filling factor This can be 
observed in Fig.|4t) for aggregates with - 0.07 (diamonds), 
where the change in filling factor shows a strong variation for 
energies above £ = lO^^A^tot^^roU- 

3.3.2. The global recipe 

In Fig.|5]a few results from the global recipe are presented, in 
which results of collisions at central impact parameter (b - 0, 
left panels) and off-center collisions (b = 0, right panels) are 
contrasted. In Figs.|5h and|5}5 the number of large particles that 
remain after a collision, A^t , is given. At low energies the number 
of fragments is the same (A^f = 1) in both cases, reflecting stick- 
ing. At very high energies {E > 5NtotEbr), the central collision 
results in catastrophic disruption (A^f - 0). Off'-center collisions, 
on the other hand, tend to produce two large fragments at higher 
energies; because they interact only with their outer parts, the 
amount of interaction is insufficient to let the colliding aggre- 
gates either stick or fragment. 

Figures |5};,d show the mass in the power-law component 
(small fragments). Central collisions result in an equal distri- 
bution of energy among the monomers. A collision energy of 
3A^tot£^br is sufficient to shatter an aggregate. Off-center colli- 
sions are more difficult to fully destroy, though, and show, more- 
over, a strong effect on porosity. In the most compact aggregate 
(crosses) over 70% of the mass ends up in the power-law compo- 
nent, whereas the remainder is in one large fragment. However, 
these are average quantities, and in some experiments all the 
mass ended up in the power-law component as can be seen from 
Fig.lSj) where A^f drops below unity. For more fluffy aggregates 
the fragmentation is much less pronounced, because the redistri- 
bution of the kinetic energy over the aggregate is less effective. 
For example, very fluffy aggregates of filling factor 0o- = 0.122 
(diamonds) colliding at an impact parameter of h = 0.75 pro- 
duce small fragments which add up to only 6% of the total mass. 
The rest of the mass is locked into two large fragments. 

The degree of damage can also be assessed through the slope 
of the power-law distribution of small fragments {q, not plot- 
ted in Fig.|5]l. The steeper the slope, the stronger the damage. 
Heavy fragmentation produces many small fragments and re- 
sults in a steepening of the power-law. Although destruction is 
very strong in the case of a central impact (the slope reaches val- 
ues of <7 = -3.7 for E > 20A^tot£^bi), it weakens considerably 



for off-center collisions (q > -2.0). For erosive events statistics 
limit an accurate determination of q. However, for erosion the 
fragments are small in any case, independent of q. 

At low energies, the amount of aggregate restructuring, 
quantified in the parameter, is independent of impact param- 
eter (Fig.|5^,f). This is simply because the collision energy is 
too low for restructuring to be significant. The aggregates' vol- 
ume then increases in a hit-and-stick fashion, resulting in a de- 
crease of the filling factor (C^ < 1). With increasing collision 
energy the degree of restructuring is enhanced, and compres- 
sion becomes more visible. Central impacts strongly affect the 
filling factor <pa-. Figure|5^ shows that the compression is max- 
imal at an impact energy of about E - A^tot^^roii- Aggregates 
that are initially compact are difficult to further compress, be- 
cause for filling factors abov e a critical value of 33% the required 
pressures increase steeply (.Blum et al.l2006l:|Paszun & Domini^ 
2008). Any further pressure will preferentially move monomers 
sideways, causing a flattening of the aggregate and a decreas- 
ing packing density. Off-center collisions, however, lead to a 
much weaker compression (Fig.|5f). Here, the forces acting on 
monomers in the impacting aggregates are more tensile, and tend 
to produce two large fragments with the filling factor unaffected, 
Q = 1. 

4. Results 

The formulation of the collision recipe in terms of the six output 
quantities enables us to calculate the collisional evolution by a 
Monte Carlo method (see Appendix|C]for its implementation). 
The sensitivity of the collisional evolution to the environment 
(e.g., gas density, grain size, grain type; see Table|2]i is assessed. 
The coagulation process is generally followed for 10^ yr. While 
we realize that bare silicates and the long timescales may not 
be fully relevant for molecular clouds, we have elected here to 
extend our calculations to fully probe the characteristics of the 
coagulation process. In particular, since fragmentation is explic- 
itly included in the collision model we evolve our runs until a 
steady-state situation materializes. 

In Sect. 14. II the results from the standard model (n = 
10^ cm"-', flo = 0.1 ijm, ice-coated silicates) are analyzed. 
Section l4~2l presents the results of our parameter study. 

4. 1 . The standard model 

Figure |6] shows the progression of the collisional evolution of 
ice-coated silicates at a density of n = 10^ cm"-' (the stan- 
dard model) starting from a monodisperse distribution of 0. 1 fim 
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Fig. 4. Quantities provided by the local recipe. The left panel shows the mass in small fragments of the power-law component, 
normalized to the reduced mass of the colliding aggregates /p„i = Mp^i/m^. The right panel shows the relative change in the 
geometrical filling factor = ^J^™/^™. Symbols refer to the initial filling factor of the larger aggregate. 




Fig. 5. Quantities provided by the global recipe. Left panels correspond to central colhsions, while the right panels correspond to 
off-center colUsions at normaUzed impact parameter 5 = 0.75. From top to bottom: number of large fragments A^f (A, B); mass of the 
small fragments component, Mpwi, normaUzed to the total mass of the two aggregates Mtot (C D); relative change in the geometrical 
filhng factor = (E, F). 
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Table 2. List of the model runs. 
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Note. (1) Model number. (2) Number density of the gas. (3) Collision 
type, describing material parameters and collision setup. Here, 'ice' 
refers to ice-coated silicates of bulk density identical to bare silicates, 
Ps = 2.65 g cm"^, but different material properties: y = 370 erg cm"^ 
and £* = 3.7 x 10'" dyn cm"^. For bare silicates, y = 25 erg cm"^ and 
£* = 2.8 X 10" dyn cm"^. (4) Monomer radius. (5) Figure reference. " 
The standard model; * filling factor of particles restricted to a minimum 
of 33%; ' central impact collisions only (b = 0). 
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Fig. 6. Mass distribution of the standard model (n 
ao = 10~^ cm, ice-coated silicates) at several times during its 
coUisional evolution, until f = 5 x 10^ yr. The distribution is 
plotted at times of 10' yr (solid lines, except for the 10^ yr curve, 
which is plotted with a dashed line) and 3 x 10' yr (all dotted 
lines), starting at f = 3 x lO'* yr. The gray shading denotes the 
spread in 10 runs. Mass is given in units of monomers. The final 
curve (thick dashed line) corresponds to 5 x 10^ yr and overlaps 
the 3 X 10^ yr curve almost everywhere, indicating that steady- 
state has been reached. 



grains. Each curve shows the average of 10 simulations, where 
the gray shading denotes the 1 cr spread in the simulations. At 
f = the distribution starts out monodisperse at size N = I. The 
distribution function /(N) gives the number of aggregates per 
unit volume such that f{N)dN is the number density of particles 
in a mass interval [A^, + dN]. Thus, at f = the initial distribu- 
tion has a number density of f(N = l,t - 0) - tijjinu/'Rgdino = 
3.5 X 10"^ cm"-' for n = 10^ cm"-' and aq = 0.1 jim. On the y-axis 
N-^f{N) is plotted, which shows the mass of the distribution per 
logarithmic interval, at several times during the collisional evo- 
lution. The mass where N^f(N) peaks is denoted the mass peak: 
it corresponds to the particles in which most of the mass is con- 
tained. The peak of the distribution curves stays on roughly the 
same level during its evolution, reflecting conservation of mass 
density. 

After 10^' yr (first solid line) a second mass peak has ap- 
peared at N ^ 10. The peak at N - 1 is a result of the compact 
(ipa- = 1) size and smaller collisional cross-section of monomers 
compared with dimers, trimers. Furthermore, the high collisional 
cross section of, e.g., dimers is somewhat overestimated, be- 
ing the result of the adopted power-law fit between the geomet- 
rical and collisional cross section (Fig. lB.3l l. These effects are 
modest, however, and do not affect the result of the subsequent 
evolution. Meanwhile, the porosity of the aggregates steadily 
increases, initially by hit-and-stick collisions but after ~10^ yr 
mostly through low-energy collisions between equal size parti- 
cles (global recipe) that do not visibly compress the aggregates. 
In Fig.|7]the porosity distribution is shown at several times dur- 
ing the collisional evolution. Initially, due to low-energy colli- 
sions the filling factor decreases as a power-law with exponent 
-0.3, 4>o- ^ A^"^--'. This trend ends after ~ 10^ at which time 
collisions have become sufficiently energetic for compaction to 
halt the fractal growth. The filling factor then stabilizes and in- 
creases only slowly. At f = 3 x 10^ yr the A^ ~ 10^ particles are 
still quite porous. 

After 3 X 10^ yr collisions have become sufficiently ener- 
getic for particles to start fragmenting, significantly changing 
the appearance of the distribution (Fig.|6]l. Slowly, particles at 
low mass are replenished and growth decelerates. When inquir- 
ing the statistics underlying the fragmenting collisions, we find 
that collisions that result in fragmentation show only modest ero- 
sion: only a tiny amount of the mass of the large aggregate is 
removed. Therefore, at the onset of erosion, growth is not imme- 
diately halted, but it is effective in replenishing the particles at 
low mass. Eventually, atN ~ 10^ (a ~ 100 fim) the erosive colli- 
sions reach a point at which there is no net growth. With increas- 
ing time and replenishment, the small particles start to reaccrete 
to produce a nearly flat distribution in terms of N-f{N). Since 
the final curve (f = 5 x 10^ yr) mostly overlaps the 3x10^ curve 
(both in Figs. |6] and |7]l it follows that steady state is reached on 
~10^ yr timescales - much longer than the timescales on which 
molecular clouds are thought to exist. 

At 10^ yr the largest particles have reached the upper limit of 
33% for the filling factor (see Fig.|7ll. Compaction increases the 
collision velocities between the particles and therefore enhances 
the fragmentation. The presence of a large population of small 
particles in the steady state distribution also hints that they are 
responsible for the higher filling factors particles of intermediate 
mass (i.e., A^ ~ 10^* - 10^) have at steady-state compared with 
the filling factor of these particles at earlier times. Indeed, the 
turnover point at A^ ~ 10^ corresponds to an energy of ^S/iioii 
these particles have with small fragments. Compaction by small 
particles is thus much more efficient than collisions with larger 
(but very fluffy) particles. 
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Fig. 7. The distribution of the filling factor, (pa-, in the standard 
model, plotted at various times. Initially, the porosity decreases 
as a power-law, ^ N^'^'^, the fractal regime. Compaction is 
most severe for the more massive particles where the filling 
factor reaches the maximum of 33%. Only mean quantities are 
shown, not the spread in (pa-. 
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time [yr] 

Fig. 8. {solid curves) The mean size (a) {dashed curves), the 
mass-weighted size {a}„, {dotted curves) and the mass-weighted 
filling factor, (0o-)m (io/Zc/ curves) of the size distribution as func- 
tion of time in the standard model {black curves), the compact 
model {dark gray curves) and the head-on only model {light gray 
curves). 



4.1.1. Compact and head-on coagulation 

To further understand the influence of the porosity on the coUi- 
sional evolution, the progression of a few key quantities as func- 
tion of time are plotted in Fig. [8] the mean size {a), the mass- 
average size {a)m, and the mass-average filhng factor {<pa-)m of 
the distribution. Here, mass-average quantities are obtained by 
weighing the particles of the Monte Carlo program by mass; e.g.. 



is the mass-weighted size. The weighing by mass has the ef- 
fect that the massive particles contribute most, because it is usu- 
ally these particles in which most of the mass resides. On the 
other hand, in a regular average all particles contribute equally, 
meaning that this quantity is particularly affected by the particles 
that dominate the number distribution. Thus, initially {a},„ = {a} 
since at f = there is only one particle size. With time, however, 
most of the mass ends up in large particles but the small particles 
still dominate by number, («),„ > (a). This picture is consistent 
with the distribution plots in Fig.|6] 

How sensitive is the emergent size distribution to the adopted 
collision recipe? To address this question we ran simulations in 
which the collision recipe is varied with respect to the standard 
model. The distribution functions of these runs are presented in 
Fig.|9] while Fig. [8] also shows the computed statistical quan- 
tities (until t - 10^ yr). In the case of compact coagulation the 
filling factor of the particles was restricted to a minimum of 33% 
(but small particles like monomers still have a higher filling fac- 
tor). Clearly, Fig. |8] shows that porous aggregates grow during 
the initial stages (cf. also Fig.|9^ and Fig.jOj^). These results are 
in line with a simple analytic model for the first stages of the 



growth, presented in Appendix |A.2| the collision timescales be- 
tween similar size aggregates is shorter when they are porous. 

Figure |9}; presents the results of the standard model in which 
collisions are restricted to take place head-on, an assumption 
that is frequently eniploye d in collision studies (e.g., Wada e t al] 
120081: ISuvama et al.ll2008l) . That is, except for the missing col- 
lision probability (/miss), the collision parameters are obtained 
exclusively from the b = entry. The temporal evolution of 
the head-on only model is also given in Fig.[8]by the light-gray 
curves. It can be seen that the particles are less porous than for 
the standard model. This follows also from the recipe, see Fig.|5] 
at intermediate energies (fi/A^tot^^roii ~ 1) central collisions are 
much more effective in compacting than off-center collisions. 
For the same reason growth in the standard model is also some- 
what faster during the early stages. However, at later times the 
differences between Fig.|9}5 and Fig.|9}; become negligible, in- 
dicating that head-on and off-center collisions do not result in a 
different fragmentation behavior. 

Thus, we conclude that inclusion of porosity significantly 
boosts the growth rates on molecular cloud relevant timescales 
(f = 10^-10^ yr). Studies that model the growth by compact par- 
ticles of the same internal density will therefore underestimate 
the aggregation. Off-center collisions are important to provide a 
(net) increase in porosity during the restructuring phase but do 
not play a critical role at later times. 

4.2. How density and material properties affect tine evolution 

Figure [TOh-c give the coUisional evolution of silicates at sev- 
eral densities. In most of the models fragmentation is impor- 
tant from the earliest timescales on. Due to the much lower 
breaking energy of silicates compared with ice, silicates already 
start fragmenting at relative velocities of ~10 m s"'. As a re- 
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Fig. 9. The effects of the colhsion recipe on the evolution of the size distribution. The standard model (b, shown for comparison) is 
varied and features: (a) compact coagulation, in which the filling factor is restricted to a lower limit of 33%; (c) head-on collisions 
only, where the impact parameter is fixed at i> = for every collision. The calculations last for 10^ yr. 



suit, the growth is very modest: only a factor of 10 in size for 
the n = 10'' cm"^ model, whereas at lower densities most of the 
mass stays in monomers. For the same reason, silicates reach 
steady state already on a timescale of 10^ yr, much faster than 
ice-coated particles. 

In the case of ice-coated silicates (Fig. [TOH-f) much higher 
energies are required to restructure and break aggregates and 
particles grow large indeed. In all cases the qualitative pic- 
ture reflects that of our standard model (Fig.fTOb). discussed in 
Sect. 14.11 porous growth in the initial stages, followed by com- 
paction and fragmentation in the form of erosion. The evolution 
towards steady-state is a rather prolonged process and is only 
complete within 10^ yr in Fig.fTOf. In the low density model 
of Fig.fTOtl fragmentation does not occur within 10^ yr. Steady 
state is characterized by a rather flat mass spectrum. 

In Fig. [TT] the collisional evolution is contrasted for three 
different monomer sizes: aq - 300 A (Fig.fTTh). 0.1 /jm (the 
standard model, Fig.fTTb). and 1 fim (Fig. [Tit). To obtain a 
proper comparison. Fig. [TT] uses physical units (grams) for the 
mass of the aggregates, rather than the dimensionless number of 
monomers, A^. From Fig.[TT]it can be seen that the models are 
extremely sensitive to the grain size. In Fig.[TTt. for example, 
the weaker aggregates result in the dominance of fragmenting 
collisions already from the start. These curves, therefore, resem- 
ble the silicate models of Fig.[TOb. 

Figure [TTh. on the other hand, shows that a reduction of the 
grain size by about a factor three (aq - 0.03 fim) enhances the 
growth significantly. Despite starting from a lower mass, the 300 
A model overtakes the standard model at f ~ lO*" yr. An under- 
standing of this behavior is provided in Appendix |A.2| the key 
element being the persistence of the hit-and-stick regime from 
which it is very difficult to break out if oq is small. Until 4x10'' yr 
visible compaction fails to take place and aggregates become 
very porous indeed (0 ^ 4 x 10""^). The consequence is that frag- 
mentation is also delayed, and has only tentatively started near 
the end of the simulations. We caution, however, against the rel- 



evance of the 300 A model; as explained in Sect. 13.1! the choice 
of flo = 300 A is too low to model aerodynamic and mechanical 
properties of MRN aggregates. But Fig. [TT]serves the purpose of 
showing the sensitivity of the collisional result on the underlying 
grain properties. 

4.3. Comparison to expected molecular cloud lifetimes 

Tables [3] and |4| present the results of the collisional evolution 
in tabular format. In Table [3] the mass- weighted size of the dis- 
tribution ((fl),„, reflecting the largest particles) is given, and in 
Table[4|the opacity of the distribution is provided, which reflects 
the behavior of the small particles. Here, opacity denotes geo- 
metrical opacity - the amount of surface area per unit mass - 
which would be applicable for visible or UV radiation, but not 
to the IR. Its definition is, accordingly. 



where the summation is over all particles in the simulation. 
These tables show, for example, that in order to grow chondrule- 
size particles (~10"^ g), dust grains need to be ice-coated and, 
except for the n = 10^ cm"-' model, coagulation times of ~10^ yr 
are needed. 

To further assess the impact of grain coagulation we must 
compare the coagulation timescales to the lifetimes of molecular 
clouds. In a study of molecular clouds in the solar neighborhood 
Hartmann et al. (2001) hint that the lifetime of molecular cloud 
is short, because of two key observational constraints: (i) most 
cores do contain young stars, rather than being starless; and (ii) 
the age of the young stars that are still embedded in a cloud is 
1-2 Myr at most. From these two arguments it follows that the 
duration of the preceding starless phase is also 1-2 Myr. If core 
lifetimes are limited to the free-fall time (Eq. ([Tji), then, the grain 
population will not leave significant imprints on either (i) the 
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Fig. 10. Distribution plots corresponding to the collisional evolution of silicates {left panels) and ice-coated silicates (right panels) 
at densities of n = 10"*, 10^ and 10^ cm"-' until t - 10^ yr. For the silicates a steady-state between coagulation and fragmentation 
is quickly established on timescales of ~10^ yr, whereas ice-coated silicates grow much larger before fragmentation kicks in. The 
initial distribution is monodisperse at oq - 10"^ cm. Note the different x-scaUng. 



large particles produced, or (ii) the removal of small particles. 
This can be seen from Tables[3]and|4]where (a),,, and (k) are also 
given at the free-fall time of the simulation (col. 6). From Table[3] 
it is seen that the sizes of the largest particles all stay below 1 fim 
(except the models that started already with a monomer sizes of 
ao - 1 jum). Likewise, Table|4]shows that the opacities from the 
tff entry are similar to those of the 'initial' 10"* yr column, i.e., 
growth is negligible on free-fall timescales. 

This information is also displayed in Fig. [12] where the opac- 
ity with respect to the initial opacity, k/kq, is plotted against 
time for all densities from the qq = 10"^ cm ice-coated sili- 



cate models. In Fig. [12] time is normalized to the initial colli- 
sion timescale between two grains, fcou.o, which is a function of 
density (see Eq. ( |A.4| i). The similarity of the curves for the first 
~10 fcou.o is in good agreement with a simple analytic model 
presented in Appendix lA.2l In models where small particles are 
replenished by fragmentation, k first obtains a minimum and 
later levels-off at k/ki) ~ 0.05. Also in Fig. [12] the free-fall and 
ambipolar diffusion timescales are indicated with diamond and 
square symbols, respectively. Due to the normalization by fcoii.o 
these occur within a relatively narrow region, despite the large 
range in densities considered. It is then clear that at a free-fall 
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Fig. 11. The effects of a different grain size ao to the collisional evolution: (a) ao = 300 yum, (b) ao - 0.1 ;um (the default, shown for 
reasons of comparison), and (c) ao = 1 A'm. To facihtate the comparison, physical units are used (grams) for the mass of aggregates, 
rather than the number of monomers (A^). 

Table 3. Mass-weighted size of the distribution, (fl),„, at several distinct events during the simulation run. 
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10^ yr 
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ftf(n) 
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1.0(-5) 


1.0(-5) 
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-5) 
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n = 




silicates 


1.0(-5) 


l.l(-5) 
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1.2(- 
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n — 


10\ 


ice 
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Note. Column (1) lists the models in terms of the density (n) and material properties. The monomer size (ao) is 0.1 fim, unless otherwise indicated. 
Cols. (2)-(5) give the mass-weighted size of the distribution at fixed coagulation times. Likewise, cols. (6)-(7) provide {a)„, at the free-fall and 
the ambipolar diffusion timescale of the cloud that corresponds to the gas density n. These are a function of density and are given in Eq. and 
Eq. respectively. Values a x 10* are denoted a(b). 



timescale no significant reduction of the opactiy takes place, 
since fff/fcoU.o ^ 1- 

However, there is still a lively debate whether the fast SF 
picture - or, rather, a short lifetime for molecular clouds - is 
general ly attainable, as cores may hav e additional support mech- 
anisms (iTassis & Mouschoviasl2004l) . If clouds are magnetically 
supported, the collapse is retarded by ambipolar diffusion (AD), 
and the relevant timescales are much longer than the free-fall 
timescale (see Eq. dU), fAD/fcoii.o s> 1 (Fig.[T2]i . Then, growth 
becomes significant, as can be seen from Table|3]where aggre- 



grates reach sizes of ~ 100 yum in the densest models on an AD- 
timescale. For the highest density models timescales are even 
sufficiently long for fragmentation to replenish the small grains. 
(Note that, although Tad decreases with density, the evolution 
of the core is determined by the quantity fAo/fcoiLO, which in- 
creases with n.) Thus, if cores evolve on AD-timescales, the ob- 
servational appearance of the core will be significantly affected. 
Tableland Fig. [T2] show that the UV-opacity, which is directly 
proportional to k, will be reduced by a factor of ~10. Studies 
that relate the Ay extinction measurements to column densities 
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Table 4. Like Table[3]but for the geometrical opacity k of the particles. 
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Fig. 12. The opacity k normalized to its initial value vs. time in 
units of the initial collision time fcoii.o (Eq. ( IA.4b ) for the ice- 
coated, ao = 0.1 ;um silicates models at five different gas densi- 
ties n. The decrease in opacity occurs on timescales of ^lOfcoii.o- 
In simulations where small grains reemerge due to fragmentation 
K starts to increase again. The free-fall {diamonds) and ambipo- 
lar diffusion timescales {squares) are indicated as far as these 
fall within 10'' yr {circles). Points of low density appear at lower 

f/fcoll,0- 



through the standard dust-to-gas ratio therefore could underesti- 
mate the amount of gas that is actually present. 



5. Discussion 

5.1. Growth characteristics and comparison to previous 
works 

In his pione ering work to the s tudy of dust coagulation in molec- 
ular clouds. Ossenkopfl (1 19931) . like our study, follows the inter- 
nal structure of particles and presents a model for the change 
in particle properties for collisions in the hit-and-stick regime. 
Furthermore, the grains are characterized by an MRN size dis- 



tribution. The model of' Ossenkopfl dl 9931) only treats the hit-and- 
stick collision regime but at the high densities (hh ^ 10* cm-^) 
and short timescales {~1Q^ yr) he considers any compaction 
or fragmentation between ice(-coated) particles is indeed of no 
concern. The coagulation then proceeds to produce particles of 
compact size ~0.5 //m at ne - 10* cm--'. It can be seen from 
Table [3]that the growth in the corresponding model of our study 
(ice, n - 10* cm-^) is higher: 2.7 jjm. This large difference 
(especially in term s of mass) can be attributed to the fact that 
FOssenkopl d 19931) ignores turbulent relative velocities between 
particles of friction times Tf < fj- As a result, growth is pre- 
dominantly PCA, because the small grains can only be swept up 
by bigger aggregates, rendering his coagulation more compact 
in comparison to our model and therefore slower. Additionally, 
due to the different definitions we use for 'density' (n vs. «h, see 
footnote |2]i our '10* cm-^ model' is denser by a factor of 1.7, 
resulting in a lower collision timescale and faster growth. 

However, at timescales f > fcoii.o (where fcoii.o for a dis- 
tribution would be the collision time between big grains) hit- 
and-stick growth will turn into CCA. Consequently, fast growth 
is expected on timescales larger than a collision timescale (see 
Appendix [All. By 10^ yr this condition has clearly been fulfilled 
in our n = 10* cm--' model, but it is likely that, due to the 
above mentioned d ifferences, i t has not been met, or perhaps 
only marginally, in lOssenkopS (Il993l) . Thus, rather than fixing 
on one point in time, a more useful comparison would be to com- 
pare the growth curves, a{t). 

On the other hand, I Weidenschilling & R uzmai kinal d 19941) . 
adopt a Bonnor-Ebert sphere to model the molecular cloud, and 
calculate the size distributi on for much lon ger tim escal e s {t = 
10^ yr). Like our study, W eidenschihing & Ruzmai kinal (1 19941) 
include fragmentation in the form of erosion and, at high ener- 
gies, shattering. Their particles are characterized by a strength 
of Q ~ 10* erg g-\ which are, therefore, somewhat weaker than 
the particles of our standard model. Although their work lacks a 
dynamic model for the porosity evolution, it is assumed that the 
initial growth follows a fractal law until 30 yum. At these sizes 
the minimum filling factor becomes less than 1%, lower than 
our results. On timescales of ~10^ yr particles grow to >100 yum, 
comparable to that of our standard model. 

A major difference between Weidensch ilhng & Ruzmaikinal 
(11994 and our works concerns the shape of the size dis- 
tribution. Whereas in our calculations the mass-peak al- 
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ways occurs at the high-mass end of the spectrum, in the 
IWeidenschilling & Ruz maikina' (1 19941) models most of the mass 
stays in th e smallest parti cles. Perhaps, the l ack of massive parti- 
cles in the IWeidenschilling & Ruzmaikinal d 19941) models is the 
result of the spatial diffusion processes this work includes; mas- 
sive particles, produced at high density, mix with less massive 
particles from the outer regions. In contrast, our findings regard- 
ing ste ady-state distributio ns agree qualitatively with the find- 
ings of lBrauer et alj ( l2008l) for protoplanetary disks. Despite the 
different environments, and therefore different velocity field, we 
find that the steady state coagulation-fragmentation mass spec- 
trum is characterized by a rather flat m^ f(m) mass function. 

It is also worthwhile to compare the aggregation results from 
our study with the constituent particles of meteorites, chon- 
drules (a ~ 300 fj.m) and calcium-aluminium inclusions (CAls, 
a ~ cm). Although most meteoriticist s accept a ne bular origin 
for these species (e.g.. lHuss et alj200ll) . lmocil(fT998 ) suggested 
that, in order to explain Al-26 free inclusions, aggregates the 
sizes of CAls (and therefore also chondrules), formed in the pro- 
tostellar cloud. These large aggregates then were self-shielded 
from the effects of the Al-26 injection event. However, our re- 
sults indicate that growth to cm-sizes seems unlikely. Only the 
dense (n > 10^) models can produce chondrule-size progenitors 
and only at a (long) ambipolar diffusion timescale. 

5.2. Observational implications for molecular clouds 

In our models we observe that the shape of the initially monodis- 
perse dust size distribution evolves first to a Gaussian-like dis- 
tribution and eventually to a flat steady-state distribution. For 
timescales longer than the coagulation timescale (Eq. (IA.4b ) we 
can expect that this result is independent of the initial condi- 
tions, even if the coagulation starts from a power-law distri- 
bution. Essentially, these distributions are a direct result of the 
physics of the coagulation: the Gaussian-like distribution reflects 
the hit-and-stick nature of the agglomeration process at low ve- 
locities while the 'flat' N^f{N) distribution at later times re- 
sults from a balance between fragmentation - erosion but not 
catastrophic destruction - and growth. In contrast, in interstel- 
lar shocks grains acquire much larger relative velocities and 
grain-grain collisions will t hen quickly shatter aggregates into 
their constituent monomers jjones et al.lll996l: iHirashita & YanI 
|2009|) . Hence, the interstellar grain size distribution will be very 
different in the dense phases of the interstellar medium than in 
the diffuse ISM and studies of the effects of grains on the opac- 
ity, ionization state and chemical inventory of molecular clouds 
will have to take this into account. 

As Fig.[T2]illustrates, in our calculations, the average surface 
area of the grain mixture - a proxy for the visual and near-IR ex- 
tinction - decreases by orders of magnitude during the initial 
coagulation process. In a general sense this finding is in agree- 
ment with observational evidence for the importance of grain 
growth in molecular clouds as obtained from studies of dust ex- 
tinction per unit column density of gas, where the latter is mea- 
sured either through HI/H2 UV absorption lines, sub-millimeter 
CO e mission lines, or X-ray absorption (cf. Whittet 2005; Jura 
Il980t IWinston et all [20071: iGoTdsmith et all 119971) . Obviously, 
this process is faster and therefore can proceed further in dense 
environments (Fig.[T2]i. As a corollary to this, the decrease in to- 
tal surface area only occurs for timescales well in excess of the 
free-fall timescale. Hence, very short lived density fluctuations 
driven by turbulences will not show this effects of coagulation 
on the total grain surface area of dust extinction. 



The study bv lChiar et alj (l2007h is - at first sight - somewhat 
at odds with this interpretation. They find that the total near-IR 
extinction keeps rising when probing deeper into dense cores 
while the strength of the 10 fim feature abruptly levels off at a 
near-IR extinction value which depends s omewhat on the cloud 
surveyed. The recent study by iMcClurd (|2009) also concludes 
that the strength of the 10 yum feature relative to the local contin- 
uum extinction decreases dramatically when the K-band extinc- 
tion exceeds 1 magnitude. Cle arly, the two featu res are carried 
by different grain populations dChiar et al.ll2007 ). Indeed, mod- 
els for interstellar extinction attribute the near-IR extinction to 
carbonaceous dust grai ns while the 10 um fea ture is a measure of 
the silicate population ( Draine & Leelll984b . Hence, these data 
suggest that silicates coagulate very rapidly when a certain den- 
sity (i.e., depth into the cloud) is reached - essentially hiding sil- 
icates grains in the densest parts of the cloud from view - while 
the carbo naceous grain population is not (as much) affected. In 
his studv. lMcClurg(l2009 ) concludes that icy grains are involved 
in this change in extinction behavior with A^. Likely, rather than 
the presence of the 13 yum ice libration band affecting the silicate 
profile, this behavior reflects grain growth. Our study shows that 
coagulation in molecular clouds is greatly assisted by the pres- 
ence of ice mantles. Once grains are covered by ice mantles, the 
increased 'stickiness' of ice takes over and the precise character- 
istics of the core become immaterial. Perhaps, therefore, the data 
suggest that silicates rapidly acquire ice mantles while carbona- 
ceous grains do not. However, there is no obvious physical basis 
for this suggestion. Further experimental studies on ice forma- 
tion on different materials will have to settle this issue. 

In this study we discussed observational implications of our 
model in a very coarse way, i.e., by considering the reduction of 
the total geometrical surface area (k) due to aggregation. We then 
find that its behavior can be largely expressed as function of the 
initial collision timescale, fcoii.o. However, for a direct compar- 
ison with observations, e.g., the 10 fim silicate absorption fea- 
ture, it is relevant to calculate the extinction properties of the 
dust distribution as function of wavelength, a nd to assess, for 
example, the significance of porous aggregates (iMin et a02008l: 
Shen et al 2008.) . This is the subject of a follow up study. 

6. Summary and conclusions 

We have studied the collisional growth and fragmentation pro- 
cess of dust in the environments of the molecular cloud (cores). 
In particular, we have focused on the collision model and the 
analysis of the several collisional evolution stages. Except for 
bouncing, the collision model features all relevant collisional 
outcomes (sticking, breakage, erosion, shattering). Furthermore, 
we have included off-center collisions in the recipe format and 
also prescribe the change to the internal structure in terms of the 
filling factor We have treated a general approach, and outcomes 
of future experiments - either numerical or laboratory - can be 
easily included. The collision model features scaling of the re- 
sults to the relevant masses and critical energies, which allows 
the calculation to proceeds beyond the sizes covered by the origi- 
nal numerical collision experiments. Our method is, in principle, 
also applicable to the dust coagulation and fragmentation stages 
in protoplanetary disks. 

We list below the key results that follow from this study: 

1 . The collisional evolution can be divided into three phases: (i) 
t < fcoU.o in which the imprints of growth are relatively mi- 
nor; (ii) a porosity-assisted growth stage, where the N^f(N) 
mass spectrum peaks at a well-defined size; and (iii) a frag- 
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mentation stage, where the N^f(N) mass spectrum is rela- 
tively flat due to the replenishment of small particles by frag- 
mentation. Fragmentation is primarily caused by erosive col- 
lisions. 

2. A large porosity speeds up the coagulation of aggregates in 
the early phases. This effect is self-enhancing, because very 
porous particles couple very well to the gas, preventing ener- 
getic collisions capable of compaction. Growth in the second 
regime can therefore become very fast. Grazing collisions 
are largely responsible for obtaining fluffy aggregates in the 
first phases, further increasing the porosity. 

3. Silicate dust grains or, in general, grains without ice-coating 
are always in the fragmentation regime. This is a result 
of their relatively low breaking energy. Freeze-out of ices, 
on the other hand, will significantly shift the fragmentation 
threshold upwards, fulfilling a prerequisite for significant ag- 
gregation in molecular clouds. 

4. Likewise, the (monodisperse) grain size that enters the col- 
lision model is also critical for the strength of the resulting 
dust aggregates. Smaller grains will increase the strength sig- 
nificantly, due to increased surface contacts. Besides, a coag- 
ulation process that starts with small grains also results in the 
creation of very porous aggregates, which further enhances 
their growth. Although a single grain size cannot fully sub- 
stitute for both the mechanical and aerodynamic properties 
of a grain size distribution, we have argued that for the MRN 
distribution a size of 0. 1 fim reflects these properties best. 

5. If cloud lifetimes are restricted to free-fall times, little co- 
agulation can be expected since the coagulation timescale 
is generally longer than However, if additional support 
mechanism are present, like ambipolar diffusion, and freeze- 
out of ice has commenced, dust aggregates of ~100 yum are 
produced, which will significantly alter the UV-opacity of 
the cloud. Conversely, our results reveal that the total dust 
surface area (and hence the extinction per H-nuclei) pro- 
vides a convenient clock that measures the lifetime of a dense 
core in terms of the initial coagulation timescale. As obser- 
vations typically reveal that the dust extinction per H-nuclei 
in dense cores has decreased substantially over that in the 
diffiise ISM, this implies that such cores are long-lived phe- 
nomena rather than transient density fluctuations. 

6. Despite the complexity of the collision model, we find that 
the decrease in (total) dust opacity can be expressed in terms 
of the initial collision time fcoU.o only, providing a relation 
for the density and lifetime of the cloud to its observational 
state (Fig. [nil. 
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Appendix A: Analytical background 

A. 1 . Relative velocities and collision timescales of dust 
particles 

The friction time, Tf , sets the amount of coupling between parti- 
cles and gas. In molecular clouds the Epstein regime is applica- 
ble for which 



Tf = 



47rCgPg cr' 



(A.l) 



where m is the mass of the particle and cr the average projected 
surface area. For compact spheres Eq. dA.ll l scales linearly with 
radius, but for porous aggregates cr can have a much steeper de- 
pendence on mass (in the case of flat structures, cr oc m) and Tf 
a much weaker dependence. For spherical grains of size oq and 
density Eq. ( lA.lb becomes 



To =Tf(flo) = 



CgPg 



(A.2) 



= 1.1 X 



-1/2 



flO 



0.1 fim j ' 



where ps = 2.65 g cm""* is used, applicable for silicates. 

Any coagulation models requires the relative velocities Au 
between two arbitrary particles. In turbulence, the motions of 
particles can become very correlated, though; e.g., particles re- 
act in similar ways to the eddy in which they are entrained. The 
mean relative motion with respect to the gas, therefore, does not 
translate to Ay. Volk et al. ( 1980) have pioneered a study to sta- 
tistically account for the collective effects of all eddies by di- 
viding the eddies into two classes - 'strong' and 'weak' - de- 
pending on the turn- over time of the eddy w ith respect to the 
particle friction time. Ormel & Cuzzi (2007) approximated the 
framework of Volketal. (1980) and Markiewicz et al. (1991) 
to provide closed-form expressions for the relative motion be- 
tween two solid particles. In general three regimes can be distin- 
guished: 

1. The low velocity regime, T2 < ti <k fj. (Here, ti > T2 are 
the friction times of the particles). Relative velocities scale 
with the absolute diff'erence in friction time. Ay oc n - T2- 

2. The intermediate velocity regime, for which fs <k ti <K Jl. 
Particle velocities scale with the square root of the largest 
particle friction time. The particle motion will not align with 
eddies of shorter turn-over time. These 'class IT eddies pro- 
vide random kicks to the particle motion - an important 
source for sustaining relative velocities of at least Au ~ u^. 

3. The heavy particle regime, T] » Jl, in which it is T2 that 
determines the relative velocity. 

Comparing the friction time of the monomer grains 
(Eq. (IA.2b ) with the smallest eddy turnover time, fj (Eq. (l6ali). 
it follows that tq > fs under most molecular cloud conditions. 



We therefore focus on the intermediate velocity regime. In par- 
ticular, the relative v elocity b etween two equal solid spheres of 
1 < To/fs < Re'^^ is dOrmel & Cuzzi..2007,) 



Avo ^ ^/3u, 



Tf 



1/2 



=8.3 X 10- cm s" 



fl() 



0.1 yum 



1/2 



(A.3) 

^-1/4/ f -.1/4 



1 105 cm-V llOK/ ■ 



Thus, velocities between silicate dust particles are ~10 m/s, and 
have a very shallow dependence on density. The same expression 
holds when the silicates are coated with ice mantles that are not 
too thick, i.e., p, is still the silicate bulk density. Dust monomers 
then collide on a collision timescale of 



Psflo'^g 



3pgA!;o 



(A.4) 



=8.5 X 10"* yr 



0.1 /im 



1/2 



ll05cm-3/ llOK/ 



where «d is the dust number density and Kgd - 100 is the stan- 
dard gas-to-dust density ratio by mass. 

Equations (IA.3I ) and ( IA.4l i are only valid for solid parti- 
cles. However, we can scale these relations to two arbitrary but 
equal aggregates of filling factor (per and (dimensionless) mass 
A^. Because m oc N and cr oc (Nl<pa-f'^ it follows that 



Tf = A^'^VrTo; 



Av^C-^ 

To 



1/2 



A^o = 



Auo; 



coll,0> 



(A.5a) 



(A.5b) 



(A.5c) 



where in the latter equation we substituted for simplicity the ge- 
ometrical cross section cr for the collisional cross-section cr'^ and 
used the monodisperse assumption oc N'^ for the dust num- 
ber density «d (this is of course only applicable for narrow dis- 
tributions). Thus, Eq. ( IA.5cb shows that the collision timescale 
decreases for very porous particles with low filling factors. 

A.2. A simple analytical model for the initial stages of the 
growth 

Despite the complexity of the recipe, it is instructive to approx- 
imate the initial collisional evolution of t he dust siz e distribu- 
tion with a simple analytical model (cf. BlumI 12004 ). Figure|7] 
suggests that the initial evolution of (pg- can be divided in two 
regimes, where the transition point occurs at a mass A^i . Initially 
(A^ < A^i), the filling factor is in the fractal regime, which can 
be well approximated by a power-law, (pg- ^ A^"-'^"'. We refer 
to the fractal regime as including hit-and-stick collisions (no re- 
structuring) as well as collisions for which E > SEyan but which 
do not lead to visible restructuring, i.e., only a small fraction of 
the grains take part in the restructuring. For N > Ni the filling 
factor starts to flatten-out. It is difficult to assign a trend for 0^ 
in the subsequent evolution. Following Fig. [T] we may assume 
that initially (per stays approximately constant for several orders 
of magnitude in A^, although at some point it will quickly as- 
sume its compact value of 33%. A sketch of the adopted porosity 
structure and the resulting scaling of velocities and timescales is 
presented in Fig. lA.ll in which it is assumed that the collapse 
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Fig. A.l. (^ray soZit/ //ne) A simplified model for the behavior of 
the filling factor with growth. Initially, for Adq < Adi, the poros- 
ity decreases (fractal growth regime). This phase is followed by 
a 'status quo' phase where filling factors will be approximately 
constant. The first compaction event is reached when velocities 
reach Ayj and fragmentation sets in when relative velocities ex- 
ceeds Ay2- {black solid line) Trend of the collision velocity and 
collision timescale. (dashed line) Trend of (Ai>)^. The numbers 
denote the power-law exponents. 



of the porous structure takes place after the point where the first 
erosive collisions occurs, atN = N2- 

From the collision recipe (Sect. lTJt we identify the critical 
energies at which visible compaction and fragmentation occur 
Compaction requires collisions between similar size particles 
(i.e., the global recipe) and Fig. |5]shows that the transition (C^ > 
1) occurs at an energy of /i/A^tot^^roU - 0.2. Regarding fragmen- 
tation, the simulations clearly show that small particles are re- 
plenished in the form of erosive collisions (local recipe). From 
Fig.|4]we assign an energy threshold of E/Nf^E\,r - 1.0. Working 
out these expressions and using a typical mass ratio of 3 for the 
global recipe (A'^ - N/6), the corresponding critical velocities 
become (Aui)^ ^ 1.0£'roll/»^o and iAv2)^ - l.QEbJmo, respec- 
tively. These energy thresholds are also indicated in Fig. lA. II 

From these expressions and the initial expressions for the rel- 
ative velocity and the collision timescale (Eqs. ( IA.3b and ( IA.4I )). 
the turn-over points A^i and N2 can be calculated. We assume 
that Avo < Avi such that a fractal growth regime exist. Then, the 
first transition point is reached at a mass 



' \Auoj \mo(Avo)^ 
. ( n 



r 



370 erg cm 



7.5 



0.1 /im 



(A.6) 



-22.5 



Unfortunately, the high powers make the numeric evaluation 
rather unstable. In our simulations we find that A^i ~ 10'*. 
Subsequently, we can write for the second transition point, the 




10 



10° 10^ 
time/fcoii,o 



10^ 



Fig. A.l. Results of the simple analytic model (solid line) and 
comparison to the {m) and {m)„, statistics of the numerical results 
of our standard model. 



onset of fragmentation, N2, 



N2 



^ . (2.0^ 

Avi / \ Eroll 



(A.7) 



: 5 X lO'' 



7 



370ergcm-2/ \0.1 //m/ \ 3.7 x 10'°dyncm 



£* 



-,-2 , 



which corresponds also well to the results from the simulation 
for which N2 ~ lOl In our simulations the first fragmentation 
involves particles that are still relatively porous, such that the 
assumption in Fig. lA. II about the porosity of the A^2-particles is 
justified. However, once steady-state has been reached, particles 
of A^2 ~ 10** will have a 33% filling factor (see Fig.|7ll. 

Using again the monodisperse assumption we also obtain 
the timescales f 1 , f2 at which these transition points are reached. 
Writing, 



dA^ 



N 

fcoll 



-1/3 



(A.8) 



where Eq. ( IA.4I ) is inserted for fcoii, we insert the power-law de- 
pendence of (pa- on A^ to obtain t. Straightforward integration 

gives 

_= A^'-'^/iSdA/'+A^-i/'O N'-^I('AN\ (A.9) 

fcollO Jl Jn, 

see Fig. lA. 21 where we plotted the results of Eq. ( |A.9t together 
with the averaged mass and the mass-averaged mass of the dis- 
tribution of the standard model. It follows that the fractal growth 
stages takes ~10 fcoii.o, or ~8 x 10^ yr (cf. ~6 x 10'' yr in the 
simulation), while A^2 is reached at ~60 fcoii.o (cf- ~30fcoii,o in the 
simulation). At larger time our fit may overestimate the growth 
rates somewhat because it assumes the filling factor stays fixed at 
its low value. Overall, the model nicely catches the trends of the 
growth and is also in good agreement with previous approaches 
(Blum 2004), although, being a monodisperse model, it cannot 
fit both (m) and (m)„,. 
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Fig.B.l. Sketch of the initial setup of our simulations. The key 
input parameters Au, b, and (pu- are illustrated. 



Appendix B: The numerical collision experiments 

The skeleton of our collision model consists of the out- 
comes of many aggregate-aggregate collision simulations. In 
this appendix we briefly review the setup of the simulations 
(Appendix lB.il ). discuss some auxiliary relations required to 
complete the collision model (Appendix lB.2l ). discuss the out- 
put format of the binary collision model (Appendix |B.3l l, and 
present a few limitations that arise due to our reliance on the 
outcomes of the numerical simulations (Appendix |B.4| ). 

B.1. Collision setup and input parameters 

Collisions between aggregates are modeled using the soft ag- 
gregates nu merical dynamics (SA ND) code (Dominik & Niibold 
I2OO2I: iPaszun & Domini3l2008h . This code treats interactions 
between indi vidual grains held to g ether by surface force s in a 
contact area dJohnson et alj Il97lt iDeriaguin et al.l [T975h . The 
SAND code calculates the equation of motion for each grain in- 
dividually and simulates vibration, rolling, twisting, and sliding 
of the grains that are in contact. These interactions lead to en- 
ergy dissipation via different channels. When two grains in con- 
tact are pulled away, the connection may break, causing loss of 
the energy. Monomers may also roll or slide over each other, 
throiigh which en ergy is also dissipated (Dominik & Tielens 
ll995L ll996lll997b . For further details regarding this model and 
testing it against laborato ry experiments we refer the reader to 
iPaszun & DominikI ( |2008|) . 

To provid e both a q ualitative and a quantitative description 
of a collision. IPaszun & Dominik ( 2009) have performed a large 
number of simulations, exploring an extensive parameter space. 
They formulate a simple collision recipe that quantifies how ki- 
netic energy, compactness, and mass ratio affect the outcome 
of aggregate-aggregate collisions. These results are presented in 
tabular format (Appendix lB.3b . Figure IbTTI sketches the setup of 
these numerical experiments, illustrating three parameters that 
shape the outcome of a collision: the initial compactness as rep- 
resented by the geometrical filling factor 0o- (see below, Eq. (|7])), 
the collision energy E, and the impact parameter b. A collision 
for each parameter set is repeated six times at different orienta- 
tions, providing information on the range of outcomes. Because 
of the occasionally fluffy structure of the aggregates not all ori- 
entations result in a collision hit, especially not those at large 
impact parameter. An overview of the parameter ranges is given 
in Table IB. Ij The radius of the monomer grains in the simula- 



Table B.l. Parameters used in the numerical simulations. 



Au 


b/b,^.^^ 


4><T 


N2IN1 


[m/s] 








(1) 


(2) 


(3) 


(4) 


0.05 


0.0 


0.070 


1.0 


0.30 


0.25 


0.090 


10-3 


0.50 


0.5 


0.122 




0.75 


0.75 


0.127 




1.0 


0.875 


0.155 




2.0 


0.95 


0.161 




4.0 




0.189 




6.0 




0.251 




8.0 








10.0 









Note. (1) relative velocity; (2) impact parameter normalized to the sum 
of the outer radii fcmax; (3) geometrical filling factor; (4) mass ratio. 



tion is flo = 0.6 jim and the adopted material properties reflect 
silicates (y = 25 erg cm"^, £ - 2.8 x 10" dyn cm"^). 

Relative velocities Av are chosen such that all relevant colli- 
sion regimes are sampled: from the pure hit-and-stick collisions, 
where particles grow without changing the internal structure of 
the colliding aggregates, up to catastrophic destruction, where 
the aggregate is shattered into small fragments. In the interme- 
diate energy regime, restructuring without fragmentation takes 
place. For aggregate collisions at large size ratios, high velocity 
impacts result in erosion of the large aggregates. 

The impact parameter b is also well sampled. We probe cen- 
tral collision {b - 0), where aggregates can be compressed, graz- 
ing impacts {b ^ ^max), where particles can be stretched due to 
inertia, and several intermediate cases. In Table lB.fj the impact 
parameter is defined relative to the outer radius of the particles, 
^max - flout, I + flout,2- Here the outer radius flout is the radius of 
the smallest sphere centered at the center-of-mass of the particle 
that fully encloses it. In the Paszun & Dominifl (12009 ) study the 
outcomes of a collision are averaged over the impact parame- 
ter b. However, in a Monte Carlo treatment, the averaging over 
impact is not necessary, and the normalized impact parameter 
b - bjbjnax can be determined using a random deviate ?, i.e., 
b - r . We have chosen to use the raw data sampled by the 
Pasz un & Do minik (2009) parameter study, explicitly including 
b as an input parameter for the Monte Carlo model. In this way 
the effects of off-center impacts can be assessed, i.e., by compar- 
ing it to models that contain only head-on collisions. 

The internal structure of the aggregates, or initial com- 
pactness, is characterized by the geometrical filling factor, 0o- 
(Eq. Q). To obtain (pa-, the projected surface area, cr, is aver- 
aged over a large number of different orientations of the particle. 
The parameter space of the filling factor (per is chosen such that 
we sample very porous , fractal ag gregates t hat grow due to the 
Brownian motion (Blum & Schrapler 20041 IPaszun & Domini^ 
I2OO6I) . through intermediate compactness aggregates that form 
by particle-cluster aggregation (PCA), up to compact aggregates 
that may result from collisional compaction. 

The final parameter that determines a collision outcome is 
the mass ratio, N2/N1 {N i being the larger aggregate). The 
IPaszun & Domini^ (|2009|) experiments sample a mass ratio be- 
tween 1 and 10 3. In this study, however, we will only use the 
collisions corresponding to the two extreme values (i.e., mass 
ratios of 1 and 10"^) as representatives of two distinct classes: 
global and local. 
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Fig. B.2. The geometrical filling factor as a function of fragment 
mass. Many simulations with different sets of parameters are 
overplotted. The dashed line indicates the least square power- 
law fit, ^ (m/moT"-^^. 



B.2. Auxiliary relations for the collision recipe 

There are a few more relations required for a consistent approach 
to the collision recipe. These are presented here for complete- 



B.2.1 . The filling factor of small fragments 

A common filling factor can be assigned to the small fragments 
produced by erosive or fragmenting collisions that constitute the 
power-law component. The compactness of these particles de- 
pends only on mass and is presented in Fig. IB. 21 where frag- 
ments produced in many simulations, reflecting a variety of col- 
lision properties, are plotted. Almost all particles are located 
along the power-law with the slope of -0.33. This provides an 
easy prescription for the filling factor of small fragments. The 
dependence indicates a fractal structure (with the fractal dimen- 
sion of about Df S3 2.0) of aggregates formed in a fragmentation 
event, since non-fractal aggregates would have a filling factor 

independent of mass^ 

As shown by Paszun & Dominil3 ( |2009|) . after reaching the 
maximum compaction, further increase of the impact energy 
produces more restructuring and results in a flattening of the 
produced aggregate. Therefore, very fluffy particles can be pro- 
duced in collisions of massive aggregates, where the power-law 
component extends to larger A^. This behavior is also observed 
in Fig. IB. 2! where fluffy, small fragments follow the power-law 
relation, while some large, still compact particles are above the 
dashed line. 



Fig. B.3. The geometrical filling factor dependence on the ratio 
of outer to geometrical radii. In this figure we plot <p,rN^-^^ to 
scale the data for aggregates of different mass. 



i.e., the collision offset is determined relative to the largest im- 
pact parameter, b,mx = flout.i +aout,2- The cross-section equivalent 
radius ao- defines our structural parameter 0,^- (see Eq. O). We 
determine the relation between the two radii (oout and ag-) em- 
pirically. Both flout and a^r are determined for many aggregates 
of various shape and mass. We sample particles with the frac- 
tal dimension in the range of Df = 1 . . . 3 and masses from 
several to a few thousands monomer masses. Th ese aggregates 
were p roduced using an algorithm developed by iFilippov et al.l 
(2000). 

Figure |B3] shows the filling factor determined for different 
aggregates versus the ratio of the outer radius over the cross- 
section equivalent radius. Diamonds correspond to the produced 
aggregates of different fractal dimension and mass. The mass de- 
pendence in Fig. |B.3| is taken into account by plotting (p,j- N^^-^^. 
In this way the data for all aggregates are well confined along a 
single curve. At small flout/flo- the curve decreases very steeply 
with increasing flout/fltr- This corresponds to compact particles 
for which flout/ao- is insensitive to filling factor. The line, how- 
ever, breaks at about flout/flo- ~ 1-2 and turns in to a power- 
law with a slope of -0.3. This shallow relation represents fluffy 
aggregates that show a large discrepancy between the projected 
surface equivalent radius and the outer radius. 

In order to provide a simple relation between the two radii, 
two power-law functions are fitted to the two regimes: compact 
particles below flout/flir = 1.2 and fluffy aggregates above that 
limit. These two functions are given by 



.compact 



^out 



0.75-4.21 log A? 



(B.la) 



B.2.2. Relation between flom and Ga- 
in this study we characterize aggregates by two different radii: 
the outer radius flout and the projected surface equivalent radius 
flo-. The first is used as a reference to the impact parameter b. 



0r^ = 1.2l(^l °'^-0.33 



(B.lb) 



To further verify these relation s we use particles produced in sev- 
eral simulations performed bv lPaszun &Dominikl(l2009h . These 
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aggregates are indicated in Fig. |B.3| bv black squares. They show 
a similar relation to the one obtained in Eq. dB.lt . Points that 
are slightly shifted above the fitted lines correspond to aggre- 
gates that are partly compressed (they did not reach the maxi- 
mum compaction). Their compact cores are still surrounded by 
a fluffy exterior that causes a small increase of the ratio of the 
outer radius over the projected surface equivalent radius Oout/flo-- 
This behavior, however, occurs at a relatively small value of 
^out/oo- < 2. At a larger size ratio the filling factor falls back 
onto the power-law given in Eq. dB.lbl ). 

We remark here that, although the fits present a general pic- 
ture, situations where aout s> a,j are not likely to material- 
ize when » 1. This would corresponds to very open frac- 
tals of fractal dimension less than two. Instead, in our mod- 
els (pa-N'^'^^ > 0.1 and we therefore always have Gqui ~ cia-- 
Consequently, the fraction of missing collisions /miss is close to 
zero in most of the cases. 

B.2.3. Hit and stick 

At very low energies {E < SZiioii) two aggregates will stick 
where they meet, without affecting the internal structure of the 
particles. This is the 'hit-and-stick' regime in which the coUi- 
sional growth can often be described by fractal laws. Two im- 
portant limits are cluster-cluster coagulation (CCA) and particle- 
cluster coagulation (PCA). In the former, two particles of equal 
size meet, which often results in very fluffy structures, whereas 
PCA describes the process in which the projectile particles are 
small with respect to the target. The filling factor then saturates 
to a constant valu e. In the case of monomers, the filling factor 
wiU reach 15% (Kozasa et al.lll992h . 

In general particles do not merely collide with either similar- 
size particles or monomers. Every size- ratio is possible and leads 
to a different change in filling factor. lOrmel et"an ( l2007h pro- 
vide an a nalytical expressio n, based upon fits to collision exper- 
iments of lOssenkopj d 19931) . that give the increase in void space 
as function of the geometrical volume of the collision partners. 
Here, the geometrical volume V is the volume that corresponds 
to the geometrical radius, Oa-. In this study additional numeri- 
cal collision experiments were used to further constrain these 
analytical fits. These experiments involved several 'monomer- 
bombardments' of aggregates of different initial filling factor. 
Using these data, we fit the volume increase as 



Table B.2. Example of an output table from the online data {fpv/i 
at = in the global recipe). 



Vyoid 

Vo 



(Vi + V2) 

N2 



36/2-1 \ 
- 1 



0.0870; 



exp 



15^2 

Vi 



0.25 



(B.2) 



where Vyoid is the increase in void space (leading to a lower (pa-), 
Vi > V2 the geometrical volumes of the collision partner, N2 the 
number of grains in the smaller particle, and Vo the monomer 
volume. The first ter m converges to CCA in the limit of V2 - Vi, 
and is the same as in lOrmel et alj (l2007h . Here, 6 - 0.95 is an 
expon ent that reflects the fractal growth in this limit dOssenkopj 
Il993h . The second expression converges to PCA in the limit of 
V2 ^ Vi. The rationale of providing a second expression is that 
in the case of V2 <K Vi (PCA) the first expression goes to zero 
very quickly (no voids are added), which is inconsistent with the 
PCA limit of 15%. From the results of our new collision exper- 
iments we have introduced an exponent of 0.25 to the PCA-part 
of Eq. ( IB. 21 ). which softens the decrease of Vvoid with increasing 
mass ratio. 







U. iZiV 


U. i J J J 


n 1 80'^ 


U.Z3UJ 


5.721(- 


-4) 


0.00000 


0.000 


0.000 


0.000 


2.595(- 


-2) 


0.00000 


2.500(-3) 


0.000 


0.000 


5.721(- 


-2) 


8.330(-4) 


0.000 


0.000 


0.000 


0.1287 




9.250(-2) 


3.042(-2) 


8.750(-3) 


2.917(-2) 


0.2288 




0.3888 


9.417(-2) 


3.042(-2) 


1.500(-2) 


0.9153 




0.9575 


0.6033 


0.2438 


0.1158 


3.661 




1.000 


1.000 


1.000 


0.7271 


8.238 




1.000 


1.000 


1.000 


1.000 


14.65 




1.000 


1.000 


1.000 


1.000 



However, not all numerical experiments could be fitted 
equally well. In fact, we had to compromise. It is likely that a 
better fit involves more parameters, e.g., the elongation of the 
aggregates or their fractal dimension. Here, we have adopted ap- 
proximate fits that follow the qualitative picture in both the CCA 
(Vi = V2) and the PCA {V2 «: Vi) limit. Remark, finally, that 
for the molecular cloud environment the hit-and-stick regime is 
only relevant in the initial stages of coagulation at densities of 
n > 10^ cm""* or grain sizes ao < 0.1 fim. 



B.3. The collision tables 

Given the level of complexity, it is not feasible to provide sim- 
ple analytical expressions for the collision outcome (in terms of 
the parameters listed in Table 13. 3T l as functio n of the collision pa- 
rameters (E, <pa-, b, Ni /A^a)- Therefore, like in lPaszun & DominikI 
(|2009), the results are expressed in a tabular format. In total 
72 tables are provided. They describe six output quantities (see 
Table l3.3] l for six impact parameters b and for both the local and 
the global recipes. Since listing all these tables here is impracti- 
cal, we will provide them in the digital form as online material 
accompanying this manuscript. We present two examples to il- 
lustrate the format. 

Each table lists one output quantity as function of the dimen- 
sionless energy parameter e and the initial filling factor of ag- 
gregates 4>tT- The only exception concerns the fraction of missed 
collisions, /Iniss- This quantity provides a correction to the col- 
lision cross-section of particles, in our case calculated from the 
outer radius of an aggregate Oout (cf. AppendixICll. The filling 
factor (per is not an appropriate quantity to use here, because it is 
ambiguous where it concerns the structure of particles. For ex- 
ample, low 4>iT could mean either a very fractal structure (and 
correspondingly high number of missing collisions) or a porous 
but homogeneous structure (and low number of missing colli- 
sions). Therefore, it is more appropriate to relate the probability 
of a collision miss to the radii with which the particle is charac- 
terized. Thus, /miss is provided as a function of the ratio of the 
outer radius over the projected surface equivalent radius, flout /Qct- 

Each table is preceded by a header that specifies: the corre- 
sponding recipe (keyword: global or local), the corresponding 
impact parameter b, and the quantity listed in the table (key- 
words are'.finiss, Nf, Sf,fpwl, q, Csig). In the case of Table|BT2] 
the header is 

# GLOBAL, b=0.0, Q=fpwl 

Therefore, Table IB .21 presents the fraction of mass in the power- 
law component, /pwi, for the global recipe and for head-on colli- 
sion. 
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Table B.3. Example of an output table from the online data (/p„i 
at = in the local recipe). 





7.009(-2) 


9.047(-2) 


0.1268 


0.1610 


0.2288 


0.000 


0.00000 


0.0000 


0.000 


0.9154 


1.001 


0.3337 


0.0000 


0.000 


3.661 


4.004 


46.55 


6.340 


0.667 


14.65 


7.007 


67.07 


35.20 


9.009 


32.95 


7.508 


148.3 


58.22 


16.02 


58.58 


9.510 


129.2 


62.40 


30.03 



In each table the first column and the first row specify the 
normalized energy parameter e and the initial filling factor 
(or the ratio of the outer over the geometrical radii flout /fl<r in the 
case of /miss), respectively. Here, s denotes the collision energy 
scaled by a normalization constant that involves (i) the breaking 
or rolling energy and (ii) the reduced or total number of par- 
ticles, see Sect. 13.21 and Table l3.3] In the case of Table [B. 2 1 the 
scaling parameter is e = E/E^rNtot- 

Table lB.3l is the second example. It is taken from the local 
recipe and it presents the /pwi quantity for the head-on collision. 
The dimensionless energy parameter s has fewer entries in the 
local recipe tables than in the global recipe. In Table |B.3l the en- 
ergy is scaled by reduced number of monomers A'^ (local recipe 
scaling) and by the breaking energy Eb^ (erosion scaling) as in- 
dicated in Table 13. 3] The header in this case is 

# LOCAL, b=0.0, Q=£pwl 

Note that in the local recipe the filling factors are lower. In this 
case larger aggregates are used to model collisions at large mass 
ratio, N1/N2 - lO-'. The fractal structure of these aggregates 
results in a lower filling factor. 

B.4. Limitations of the collision recipe 

The main limitation of the collision recipe is that, due to com- 
putational constraints, the binary aggregate simulations can only 
simulate aggregates of a mass < lO-'. For the recipe to be- 
come applicable for large aggregates scaling of the results of the 
collision experiments is required. This is a critical point of the 
recipe for which suitable dimensionless quantities had to be de- 
termined. However, the extrapolation assumes that the collision 
physics that determines the outcomes of collisions at low-A^ also 
holds at large scales. This is a crucial assumption in which colli- 
sional outcomes like bouncing are a priori not possible because 
these do not take place at the low-A^ part of the simulations. 

Bou ncing of aggreg ates is observed in laboratory exper- 
iments (Blum & MUncW 'l993t iBlumI 120061: iLangkows ki et all 
[2OO8; Weidling et al. 2009), whereas it does not occur in our 
simulations. For silicates, bouncing occurs at sizes above ap- 
proximately 100 fim (i.e., > lO** particles) and is not fully 
understood from a microphysical perspective. In the case of ice- 
coated silicate grains, which provide stronger adhesion forces, 
our simulations show that growth proceeds to ~100 fim sizes. 
In this case, therefore, bouncing might slow down the growth 
earlier than observed in our experiments, especially when the 
internal structure has aheady re-adjusted to a compact state. 
However, it is presently unclear how these laboratory experi- 
ments apply to ice aggregates and hence whether and to what 
extent the results would be affected by bouncing. We recognize 
that this may, potentially, present a limitation to the growth of 



aggregates in molecular clouds, but also emphasize it will not af- 
fect the main conclusions from this study as in only a few models 
aggregates grow to sizes »100 //m. 

Another assumption of the collision model is that the grains 
have a spherical geometry. Again, computational constraints rule 
out numerical modeling of randomly shaped particles. Whether 
erratically shaped grains would help or harm the sticking or 
bouncing is unclear Because the strength of an aggregate is de- 
termined by the amount of contact area between two grains, the 
strength of irregularly shaped monomers depends on the local 
radius of curvature. Therefore, highly irregular grains are held 
by contacts of much smaller size, because they are connected 
by surface asperities. On the other hand, irregular grains may 
form more than one contact. However, the geometry of the grains 
does not necessarily pose a bottleneck to the validity of the col- 
lision model. Instead, like the size distribution, the consequence 
of irregularly shaped monomers is reflected in a different energy 
scaling. 

Appendix C: The Monte Carlo program 

The advantage of a Monte Carlo (MC) approach to the calcula- 
tion of the collisional evolution is that collisions are modeled 
individually and that they, therefore, bear a direct correspon- 
dence to the collision model. Furthermore, in a MC approach 
structural parameters (like (pa-) can be easily included and the 
collisional outcome can be quantified in detail. From the two 
particle properties (A^ and 0o-) the collisional quantities are de- 
rived, e.g., the relative velocities Av between the aggregates (see 
Appendix |A. It . Then, from Au and the particle's outer radii we 
calculate the collision rates between all particles present in the 
MC simulation. After the MC model has selected the collision 
partners, the collision recipe is implemented. First, the particle 
properties {m,(pa-) and the collision properties {Au) are turned 
into a collision 'grid point' given by the dimensionless e, ^o- and 
b. The six collision quantities (Table [33) are then taken from the 
appropriate entries from the recipe tables. Finally, these quanti- 
ties specify the change to the initial particle properties (m, (pa-) 
and also describe the properties of the collision fragments. 

By virtue of the scaling relations discussed in Sect. |3.2| the 
MC coagulation model is able to treat much larger aggregates 
than the binary collision experiments. A broad size distribu- 
tions, which may, e.g., result from injection of small particles 
due to fragmentation, can, however, become problematic for a 
MC approach, since the high dynamic range required consumes 
computational resources. To o vercome this prob l em w e use the 
grouping method outlined by lOrmel & SpaansI (l2008l) . In this 
method the 1-1 correspondence between a simulation particle 
and a physical particle is dropped; instead, the simulated par- 
ticles are represented by groups of identical physical particles. 
The group's mass is determined by the peak of the m^ f(m) mass 
distribution - denoted nip - and particles of smaller mass 'travel' 
together in groups of total mass nip. Grouping entails that a large 
particle can collide with many small particles simultaneously - 
a necessary approximation of the collision process. 

Below, we present the way in which we have imple- 
mented the collision recipe with the grouping method of 
IOrmel& SpaansI (I2OO8I) . 

C.1. Collision rates 

The cycle starts with the calculation (or update) of the collision 
rates between the groups of the simulation. The individual col- 
lision rate between two particles / and j is dj - KijfV (units: 
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s '), where is the simulation volume and K the collision ker- 
nel. For grouped collisions C,y is larger because many particles 
are involved in the collision. The collision kernel K is defined 
as Kij - (/rAvii with cr?. = :'!'(aouti + «out2)^ the collisional 
cross section (uncorrected for missing collisions) and Ai;,j the 
average root-mean-square relative velocity (See Appendix lA.il ). 
Thus, to calculate the collision rates we need the relative veloci- 
ties and the relation between the geometrical and the outer radius 
(Appendix |B.2.2| i. 

C.2. Determination of collision partners 

Random numbers determine which two groups collide and the 
number of particles that are involved from the / and /' groups, rji 
and rjj. Then, each /-particle collides with rjj/rji j-particles. The 
grouping method implicitly assumes that collision rates do not 
change significantly during the collision process. To enforce the 
plausibility of this assumption the grouping method limits the 
total mass of the j-particles colliding with the /-particle to be at 
most 1% of the mass of an /-particle, i.e., rijiiijlrjimi < fc - 10"^. 
Therefore, grouped collisions occur only in the local recipe. For 
erosion or sticking this procedure works as intended. However, 
in collisions that result in breakage the grouping assumption is 
potentially problematic, since the particle properties - and hence 
the collision rates - then clearly change significantly over a sin- 
gle collision. Fortunately, in the local recipe breakage is rela- 
tively unimportant. Catastrophic disruptions (shattering) is prob- 
lematic for the same reasons, because when it occurs, there is no 
'large' aggregate left. However, for energetic reasons we expect 
that shattering occurs mainly when two equal size particles are 
involved, for which the global recipe would apply (and no group- 
ing). In the following we continue with a collision of t], - rjj/rji 
7-particles colliding with a single /-particle. 

C.3. Determining the collision quantities in grouped collisions 

When the collision is in the 'hit-and-stick' regime the properties 
of the new particles are easily found by adding the masses of 
the J-particles to the /-particle and calculating their filling fac- 
tor using Eq. ( IB. 2b . We therefore concentrate here on the local 
or global recipe. The collision is then characterized by the three 
dimensionless parameters: normalized collision energy e, filling 
factor 0O- and impact parameters b (Sect. [J!2T i. These three pa- 
rameters constitute an arbitrary point in the 3D {s, (pg-, b)-space, 
and will in general be confined by eight grid points (k) which 
correspond to the parameters at which results from the collision 
experiments are available, see Fig. lC.ll We next distribute the 
rjt collisions over the grid point in which the weight of a grid 
point is inversely proportional to the 'distance' to (e, <p,j, b) as 
explained in Fig. lC.ll Taking account of the collisions that re- 
sult in a miss, we have 

8 8 

nt - '7miss + ^ ^7*; '?miss - ^ '7miss,A:, (C.l) 

k=\ k=\ 

where T/miss.* - ?/;^'A/miss,t denotes the number of collisions at 
the grid point resulting in a miss. Here, Pk denotes the weight of 
the grid point i^^. Pk - 1), /miss the fraction of missed collisions 
at the grid point and the ^ sign indicates this number is rounded 
to integer values. Similarly, the number of 'hits' at a grid point is 
given by ^ rjiPki^ - /miss,*-)- Not all of these grid points have 
to be occupied (i.e., rjit can be zero). In the special case without 
grouping rj, = I and one grid point at most is occupied. 




£ - £1 £2 - £ 



Fig. C.l. Illustration of the picking of the grid points. The colli- 
sion takes place at (e, (pa-, b): a point that is generally surrounded 
by eight grid points (here corresponding to the nodes of the cube) 
at which results from the binary collision simulations are avail- 
able. Each node is then assigned a probability inversely propor- 
tional to the distance to the grid point. Thus, the probability that 
the energy parameter e = si is picked (corresponding to four of 
the eight grid nodes) is Pi = (e - ei)/(e2 - si)- The procedure is 
identical for the other grid points. 

We continue with the general case of multiple occupied grid 
points. First, we consider the mass that is eroded, given by the 
fpv/\,k quantities. The mass eroded at one grid point per collision 
is given by Mpwi,;t = U-^ikim + mj) (global recipe) or Mpwi.t = 
f^viik'nimjlinii + mj) (local recipe). Then, the total mass eroded 
by the group collision is 

8 

Mpwi = ^ Mpwi,^?/^ (C.2) 

If this quantity is more than m,, clearly there is no large fragment 
component|3 Otherwise, the mass of the large fragment compo- 
nent is Miarge = '«/ + {rj, - ri^i^^)m j - Mpwi. Each Mpwi,* quan- 
tity is distributed as a power-law with the exponent provided by 
the slope of the grid point (see below). Concerning the large- 
fragment component, there is a probability that it will break, 
given by the Nf^k andSf t quantities. As argued before, breakage 
within the context of the grouping algorithm cannot be consis- 
tently modeled. Notwithstanding these concerns, we choose to 
implement it in the grouping method. Because its probability is 
small, we assume it happens at most only once during the group 
collision. The probability that it occurs is then 

8 

P2^l-Y](^-P2,k)"', (C.3) 

k=l 

where P2,a is the probability that breakage occurs at a grid point 
and follows from the Sf and A^f quantities. If breakage occurs, 

' Recall that in grouped collisions (rjt > 1) this implies that the group- 
ing method is not fully accurate as the change in mass is of the order of 
the mass itself; but the procedure is always fine if rj, = I. 
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the masses Mp„i are removed first and we divide the remaining 
mass Miai-ge in two. 

The last quantity to determine is the change in the filling 
factor of the large aggregate, denoted by the symbol for one 
collision. Like Eq. ( |C.3t we multiply the changes in at the 
individual grid nodes, 

8 

0<r,large = <^*<r)«, , (C.4) 

k=\ 

This completes the implementation of the coUisional outcome 
within the framework of the grouping mechanism. That is, we 
have the masses and the filling factor of the large fragment com- 
ponent (Miarge, </i(7-,iai-ge), and havc computed the distribution of 
the power-law component in terms of mass. Recall, finally, that 
all these results are per /-particle, and that the multiplicity of the 
results is 77,. 



I2008h to combine neighboring species, a process that averages 
over their (structural) parameters but conserves mass. This sig- 
nificantly improved the efficiency (i.e., speed) of the simulation, 
although the many fragments created by the collisions (all con- 
tributing to a higher A^s) can be regarded as a redundancy, be- 
cause it requires a lot of subsequent regrouping. The alterna- 
tive would be to produce only 1 new species per collision event 
(see Zsom & Dullemond 2008). Here, we prefer to stick with a 
more detailed representation of each collision event by creating a 
range of particles, but we acknowledge that this amount of detail 
is to some extent lost by the subsequent merging. 



C.4. Picking of ttie power-law component masses 

The final part of the MC cycle is to pick particles according to 
the power-law distribution, under the constraints of a total mass 
IkMpviik and slope qi^ at each grid point k. The general formula 
for picking the mass of the small fragments is 



'"small = mo 



1 +? 



^rem 



mo 



l+q \ 
- 1 



1/(1+?) 

(C.5) 



where mi-em is the remaining mass of the distribution (it starts 
at mrem = T]kMp„iit and decreases every step by msmaii) and ? 
a random number between and 1 . From the definition of the 
power-law component msmaii cannot be more than 25% of the 
total mass. In the MC program the number of distinct fragments 
that can be produced is limited to a few per grid point. This is to 
prevent an influx of a very large number of species (non-identical 
particles; in this case, particles of different mass), which would 
lead to severe computational problems, filling-up the statevector 
array (see below). Therefore, if the same mass msmaii is picked 
again it is considered to be the same species, and the multiplicity 
of this species is increased by one. After we have obtained a 
maximum of //dis distinct species, we redistribute the mass Mpow 
over the species. In this way the fragment distribution is only 
sampled at a few discrete points. 



C.5. Merging/Duplication 

The final part of the MC program consist of an inventory, and 
possible adjustment, of t he amount of groups an d species (A^s) 
present in the program jOrmel & SpaansI 120081) . To combine 
a sufficiently high resolution with an efficient computation in 
terms of speed is one of the virtues of the grouping method. 
One key parameter, determining the resolution of the simula- 
tion, is the A^* parameter (the target number of species in a sim- 
ulation). In order to obtain a sufficient resolution we require that 
a total mass of mp(t)N* is present in the simulation at all times, 
where mp(r) is the mass peak of the distribution, nij, = {m^}/{m}. 
Particles are duplicated to fulfill this criterion, adding mass to 
the system. To prevent a pileup of spec ies we followed the 'equal 
mass method ' ilOrmel & Spamsll2008l) . However, we found that 
due to the fragmentation many species were created at any 
rate - too many, in fact (A^s > N*) which would severely af- 
fect the efficiency of the program. Therefore, when A^s - 2N* 
was reached we used the 'merging algorithm' dOrmel & SpaansI 
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Appendix D: List of symbols 



Symbol Description 





Reduced modulus of elasticity (Eq. (|9Jl) 




Gas-to-dust ratio by mass 


Au 


Relative velocity (AppendixlAt 


7 


Surface energy density (Eq. (|9]l) 




Number of particles or groups (AppendixICli 


(pa 


Geometrical filling factor (Eq. (|7]i) 




Molecular mass (Sect.|2]i 


Vm,Vt 


Molecular, turbulent viscosity (Sect. 12.11 


fcrit 


Critical displacement for irreversible rolling (Eq. ^) 


Ps 


Material density, ps = 2.65 g cm""* (silicates) 


Pg 


Gas density, pg - junniu 


cr 


Average projected surface area (Sect.|3ll 




Collisional cross section (Sect.|3]l 


Tf 


Friction time (Eq. (lA.lb) 




Change in geometrical filling factor, = <pcr/<p'^^ 




(Sect.[33]l 


Df 


Fractal dimension 


E 


Collision energy, E = im^(Ay)^ 


Emll 


Rolling energy (Eq. (|9|) 


Ehr 


Breaking energy (Eq. (|9])) 


N 


Number of grains in aggregate (dimensionless measure 




of mass) 




Reduced number of monomers in collision A'^ = 




NiN2/iNi +N2) 


Ni 


Number of big fragments 


Mot 


Total number of monomers in collision, A^tot = A^i + N2 


Re 


Reynolds number (Eq. Q) 


Sf 


Spread in number of fragments of big component 




(Sect. [ID 


St 


Particle Stokes number (AppendixlAb 


T 


Temperature (Sect. 12.11) 


flO 


Monomer radius 


^out 


Aggregate outer radius (Fig.[T]i 


flo- 


Aggregate geometrical radius (projected surface equiva- 




lent radius) (Fig.lTJ 




Reduced radius (Eq. (|9]l) 


b 


Impact parameter 


^max 


Sum of the outer radii, b^^^ - ai.out + fla.out 


b 


Normalized impact parameter, b = b/b^^x 


Cg 


Sound speed (gas) 


/miss 


Fraction of collision misses (Sect. 13.3b 


/pwl 


Fraction of mass in power-law component (Sect. [33) 


« 


Particle density (gas) 


m 


Particle mass 




Reduced mass 


mil 


Hydrogen mass 


q 


Power-law exponent (size distribution) (Sect. 13. 3b 


r 


Random deviate 


hi 


Ambipolar diffusion time (Eq. ([2]l) 


fcoU.O 


Initial collision time (Eq. (|A.4b) 


?fr 


Free-fall time (Eq. ([B) 


U 


Inner (Kolmogorov) eddy turn-over time (Eq. (l6ab) 


ul 


Large eddy turn-over velocity (Sect. 12. lb 




Inner (Kolmogorov) eddy turn-over velocity (Eq. (|6bli) 



